source: trunk/Clp/src/ClpConstraintQuadratic.cpp @ 1525

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

  • Property svn:keywords set to Id
File size: 10.3 KB
Line 
1/* $Id: ClpConstraintQuadratic.cpp 1525 2010-02-26 17:27:59Z mjs $ */
2// Copyright (C) 2007, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5#include "CoinPragma.hpp"
6#include "CoinHelperFunctions.hpp"
7#include "CoinIndexedVector.hpp"
8#include "ClpSimplex.hpp"
9#include "ClpConstraintQuadratic.hpp"
10#include "CoinSort.hpp"
11//#############################################################################
12// Constructors / Destructor / Assignment
13//#############################################################################
14
15//-------------------------------------------------------------------
16// Default Constructor
17//-------------------------------------------------------------------
18ClpConstraintQuadratic::ClpConstraintQuadratic ()
19     : ClpConstraint()
20{
21     type_ = 0;
22     start_ = NULL;
23     column_ = NULL;
24     coefficient_ = NULL;
25     numberColumns_ = 0;
26     numberCoefficients_ = 0;
27     numberQuadraticColumns_ = 0;
28}
29
30//-------------------------------------------------------------------
31// Useful Constructor
32//-------------------------------------------------------------------
33ClpConstraintQuadratic::ClpConstraintQuadratic (int row, int numberQuadraticColumns ,
34          int numberColumns, const CoinBigIndex * start,
35          const int * column, const double * coefficient)
36     : ClpConstraint()
37{
38     type_ = 0;
39     rowNumber_ = row;
40     numberColumns_ = numberColumns;
41     numberQuadraticColumns_ = numberQuadraticColumns;
42     start_ = CoinCopyOfArray(start, numberQuadraticColumns + 1);
43     int numberElements = start_[numberQuadraticColumns_];
44     column_ = CoinCopyOfArray(column, numberElements);
45     coefficient_ = CoinCopyOfArray(coefficient, numberElements);
46     char * mark = new char [numberQuadraticColumns_];
47     memset(mark, 0, numberQuadraticColumns_);
48     int iColumn;
49     for (iColumn = 0; iColumn < numberQuadraticColumns_; iColumn++) {
50          CoinBigIndex j;
51          for (j = start_[iColumn]; j < start_[iColumn+1]; j++) {
52               int jColumn = column_[j];
53               if (jColumn >= 0) {
54                    assert (jColumn < numberQuadraticColumns_);
55                    mark[jColumn] = 1;
56               }
57               mark[iColumn] = 1;
58          }
59     }
60     numberCoefficients_ = 0;
61     for (iColumn = 0; iColumn < numberQuadraticColumns_; iColumn++) {
62          if (mark[iColumn])
63               numberCoefficients_++;
64     }
65     delete [] mark;
66}
67
68//-------------------------------------------------------------------
69// Copy constructor
70//-------------------------------------------------------------------
71ClpConstraintQuadratic::ClpConstraintQuadratic (const ClpConstraintQuadratic & rhs)
72     : ClpConstraint(rhs)
73{
74     numberColumns_ = rhs.numberColumns_;
75     numberCoefficients_ = rhs.numberCoefficients_;
76     numberQuadraticColumns_ = rhs.numberQuadraticColumns_;
77     start_ = CoinCopyOfArray(rhs.start_, numberQuadraticColumns_ + 1);
78     int numberElements = start_[numberQuadraticColumns_];
79     column_ = CoinCopyOfArray(rhs.column_, numberElements);
80     coefficient_ = CoinCopyOfArray(rhs.coefficient_, numberElements);
81}
82
83
84//-------------------------------------------------------------------
85// Destructor
86//-------------------------------------------------------------------
87ClpConstraintQuadratic::~ClpConstraintQuadratic ()
88{
89     delete [] start_;
90     delete [] column_;
91     delete [] coefficient_;
92}
93
94//----------------------------------------------------------------
95// Assignment operator
96//-------------------------------------------------------------------
97ClpConstraintQuadratic &
98ClpConstraintQuadratic::operator=(const ClpConstraintQuadratic& rhs)
99{
100     if (this != &rhs) {
101          delete [] start_;
102          delete [] column_;
103          delete [] coefficient_;
104          numberColumns_ = rhs.numberColumns_;
105          numberCoefficients_ = rhs.numberCoefficients_;
106          numberQuadraticColumns_ = rhs.numberQuadraticColumns_;
107          start_ = CoinCopyOfArray(rhs.start_, numberQuadraticColumns_ + 1);
108          int numberElements = start_[numberQuadraticColumns_];
109          column_ = CoinCopyOfArray(rhs.column_, numberElements);
110          coefficient_ = CoinCopyOfArray(rhs.coefficient_, numberElements);
111     }
112     return *this;
113}
114//-------------------------------------------------------------------
115// Clone
116//-------------------------------------------------------------------
117ClpConstraint * ClpConstraintQuadratic::clone() const
118{
119     return new ClpConstraintQuadratic(*this);
120}
121
122// Returns gradient
123int
124ClpConstraintQuadratic::gradient(const ClpSimplex * model,
125                                 const double * solution,
126                                 double * gradient,
127                                 double & functionValue,
128                                 double & offset,
129                                 bool useScaling,
130                                 bool refresh) const
131{
132     if (refresh || !lastGradient_) {
133          offset_ = 0.0;
134          functionValue_ = 0.0;
135          if (!lastGradient_)
136               lastGradient_ = new double[numberColumns_];
137          CoinZeroN(lastGradient_, numberColumns_);
138          bool scaling = (model && model->rowScale() && useScaling);
139          if (!scaling) {
140               int iColumn;
141               for (iColumn = 0; iColumn < numberQuadraticColumns_; iColumn++) {
142                    double valueI = solution[iColumn];
143                    CoinBigIndex j;
144                    for (j = start_[iColumn]; j < start_[iColumn+1]; j++) {
145                         int jColumn = column_[j];
146                         if (jColumn >= 0) {
147                              double valueJ = solution[jColumn];
148                              double elementValue = coefficient_[j];
149                              if (iColumn != jColumn) {
150                                   offset_ -= valueI * valueJ * elementValue;
151                                   double gradientI = valueJ * elementValue;
152                                   double gradientJ = valueI * elementValue;
153                                   lastGradient_[iColumn] += gradientI;
154                                   lastGradient_[jColumn] += gradientJ;
155                              } else {
156                                   offset_ -= 0.5 * valueI * valueI * elementValue;
157                                   double gradientI = valueI * elementValue;
158                                   lastGradient_[iColumn] += gradientI;
159                              }
160                         } else {
161                              // linear part
162                              lastGradient_[iColumn] += coefficient_[j];
163                              functionValue_ += valueI * coefficient_[j];
164                         }
165                    }
166               }
167               functionValue_ -= offset_;
168          } else {
169               abort();
170               // do scaling
171               const double * columnScale = model->columnScale();
172               for (int i = 0; i < numberCoefficients_; i++) {
173                    int iColumn = column_[i];
174                    double value = solution[iColumn]; // already scaled
175                    double coefficient = coefficient_[i] * columnScale[iColumn];
176                    functionValue_ += value * coefficient;
177                    lastGradient_[iColumn] = coefficient;
178               }
179          }
180     }
181     functionValue = functionValue_;
182     offset = offset_;
183     CoinMemcpyN(lastGradient_, numberColumns_, gradient);
184     return 0;
185}
186// Resize constraint
187void
188ClpConstraintQuadratic::resize(int newNumberColumns)
189{
190     if (numberColumns_ != newNumberColumns) {
191          abort();
192#ifndef NDEBUG
193          int lastColumn = column_[numberCoefficients_-1];
194#endif
195          assert (newNumberColumns > lastColumn);
196          delete [] lastGradient_;
197          lastGradient_ = NULL;
198          numberColumns_ = newNumberColumns;
199     }
200}
201// Delete columns in  constraint
202void
203ClpConstraintQuadratic::deleteSome(int numberToDelete, const int * which)
204{
205     if (numberToDelete) {
206          abort();
207          int i ;
208          char * deleted = new char[numberColumns_];
209          memset(deleted, 0, numberColumns_ * sizeof(char));
210          for (i = 0; i < numberToDelete; i++) {
211               int j = which[i];
212               if (j >= 0 && j < numberColumns_ && !deleted[j]) {
213                    deleted[j] = 1;
214               }
215          }
216          int n = 0;
217          for (i = 0; i < numberCoefficients_; i++) {
218               int iColumn = column_[i];
219               if (!deleted[iColumn]) {
220                    column_[n] = iColumn;
221                    coefficient_[n++] = coefficient_[i];
222               }
223          }
224          numberCoefficients_ = n;
225     }
226}
227// Scale constraint
228void
229ClpConstraintQuadratic::reallyScale(const double * )
230{
231     abort();
232}
233/* Given a zeroed array sets nonquadratic columns to 1.
234   Returns number of nonlinear columns
235*/
236int
237ClpConstraintQuadratic::markNonlinear(char * which) const
238{
239     int iColumn;
240     for (iColumn = 0; iColumn < numberQuadraticColumns_; iColumn++) {
241          CoinBigIndex j;
242          for (j = start_[iColumn]; j < start_[iColumn+1]; j++) {
243               int jColumn = column_[j];
244               if (jColumn >= 0) {
245                    assert (jColumn < numberQuadraticColumns_);
246                    which[jColumn] = 1;
247                    which[iColumn] = 1;
248               }
249          }
250     }
251     int numberCoefficients = 0;
252     for (iColumn = 0; iColumn < numberQuadraticColumns_; iColumn++) {
253          if (which[iColumn])
254               numberCoefficients++;
255     }
256     return numberCoefficients;
257}
258/* Given a zeroed array sets possible nonzero coefficients to 1.
259   Returns number of nonzeros
260*/
261int
262ClpConstraintQuadratic::markNonzero(char * which) const
263{
264     int iColumn;
265     for (iColumn = 0; iColumn < numberQuadraticColumns_; iColumn++) {
266          CoinBigIndex j;
267          for (j = start_[iColumn]; j < start_[iColumn+1]; j++) {
268               int jColumn = column_[j];
269               if (jColumn >= 0) {
270                    assert (jColumn < numberQuadraticColumns_);
271                    which[jColumn] = 1;
272               }
273               which[iColumn] = 1;
274          }
275     }
276     int numberCoefficients = 0;
277     for (iColumn = 0; iColumn < numberQuadraticColumns_; iColumn++) {
278          if (which[iColumn])
279               numberCoefficients++;
280     }
281     return numberCoefficients;
282}
283// Number of coefficients
284int
285ClpConstraintQuadratic::numberCoefficients() const
286{
287     return numberCoefficients_;
288}
Note: See TracBrowser for help on using the repository browser.