Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (4 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpQuadraticObjective.cpp

    r1723 r2385  
    1616// Default Constructor
    1717//-------------------------------------------------------------------
    18 ClpQuadraticObjective::ClpQuadraticObjective ()
    19      : ClpObjective()
    20 {
    21      type_ = 2;
    22      objective_ = NULL;
    23      quadraticObjective_ = NULL;
    24      gradient_ = NULL;
    25      numberColumns_ = 0;
    26      numberExtendedColumns_ = 0;
    27      activated_ = 0;
    28      fullMatrix_ = false;
     18ClpQuadraticObjective::ClpQuadraticObjective()
     19  : ClpObjective()
     20{
     21  type_ = 2;
     22  objective_ = NULL;
     23  quadraticObjective_ = NULL;
     24  gradient_ = NULL;
     25  numberColumns_ = 0;
     26  numberExtendedColumns_ = 0;
     27  activated_ = 0;
     28  fullMatrix_ = false;
    2929}
    3030
     
    3232// Useful Constructor
    3333//-------------------------------------------------------------------
    34 ClpQuadraticObjective::ClpQuadraticObjective (const double * objective ,
    35           int numberColumns,
    36           const CoinBigIndex * start,
    37           const int * column, const double * element,
    38           int numberExtendedColumns)
    39      : ClpObjective()
    40 {
    41      type_ = 2;
    42      numberColumns_ = numberColumns;
    43      if (numberExtendedColumns >= 0)
    44           numberExtendedColumns_ = CoinMax(numberColumns_, numberExtendedColumns);
    45      else
    46           numberExtendedColumns_ = numberColumns_;
    47      if (objective) {
    48           objective_ = new double [numberExtendedColumns_];
    49           CoinMemcpyN(objective, numberColumns_, objective_);
    50           memset(objective_ + numberColumns_, 0, (numberExtendedColumns_ - numberColumns_)*sizeof(double));
    51      } else {
    52           objective_ = new double [numberExtendedColumns_];
    53           memset(objective_, 0, numberExtendedColumns_ * sizeof(double));
    54      }
    55      if (start)
    56           quadraticObjective_ = new CoinPackedMatrix(true, numberColumns, numberColumns,
    57                     start[numberColumns], element, column, start, NULL);
    58      else
    59           quadraticObjective_ = NULL;
    60      gradient_ = NULL;
    61      activated_ = 1;
    62      fullMatrix_ = false;
     34ClpQuadraticObjective::ClpQuadraticObjective(const double *objective,
     35  int numberColumns,
     36  const CoinBigIndex *start,
     37  const int *column, const double *element,
     38  int numberExtendedColumns)
     39  : ClpObjective()
     40{
     41  type_ = 2;
     42  numberColumns_ = numberColumns;
     43  if (numberExtendedColumns >= 0)
     44    numberExtendedColumns_ = CoinMax(numberColumns_, numberExtendedColumns);
     45  else
     46    numberExtendedColumns_ = numberColumns_;
     47  if (objective) {
     48    objective_ = new double[numberExtendedColumns_];
     49    CoinMemcpyN(objective, numberColumns_, objective_);
     50    memset(objective_ + numberColumns_, 0, (numberExtendedColumns_ - numberColumns_) * sizeof(double));
     51  } else {
     52    objective_ = new double[numberExtendedColumns_];
     53    memset(objective_, 0, numberExtendedColumns_ * sizeof(double));
     54  }
     55  if (start)
     56    quadraticObjective_ = new CoinPackedMatrix(true, numberColumns, numberColumns,
     57      start[numberColumns], element, column, start, NULL);
     58  else
     59    quadraticObjective_ = NULL;
     60  gradient_ = NULL;
     61  activated_ = 1;
     62  fullMatrix_ = false;
    6363}
    6464
     
    6666// Copy constructor
    6767//-------------------------------------------------------------------
    68 ClpQuadraticObjective::ClpQuadraticObjective (const ClpQuadraticObjective & rhs,
    69           int type)
    70      : ClpObjective(rhs)
    71 {
    72      numberColumns_ = rhs.numberColumns_;
    73      numberExtendedColumns_ = rhs.numberExtendedColumns_;
    74      fullMatrix_ = rhs.fullMatrix_;
    75      if (rhs.objective_) {
    76           objective_ = new double [numberExtendedColumns_];
    77           CoinMemcpyN(rhs.objective_, numberExtendedColumns_, objective_);
    78      } else {
    79           objective_ = NULL;
    80      }
    81      if (rhs.gradient_) {
    82           gradient_ = new double [numberExtendedColumns_];
    83           CoinMemcpyN(rhs.gradient_, numberExtendedColumns_, gradient_);
    84      } else {
    85           gradient_ = NULL;
    86      }
    87      if (rhs.quadraticObjective_) {
    88           // see what type of matrix wanted
    89           if (type == 0) {
    90                // just copy
    91                quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
    92           } else if (type == 1) {
    93                // expand to full symmetric
    94                fullMatrix_ = true;
    95                const int * columnQuadratic1 = rhs.quadraticObjective_->getIndices();
    96                const CoinBigIndex * columnQuadraticStart1 = rhs.quadraticObjective_->getVectorStarts();
    97                const int * columnQuadraticLength1 = rhs.quadraticObjective_->getVectorLengths();
    98                const double * quadraticElement1 = rhs.quadraticObjective_->getElements();
    99                CoinBigIndex * columnQuadraticStart2 = new CoinBigIndex [numberExtendedColumns_+1];
    100                int * columnQuadraticLength2 = new int [numberExtendedColumns_];
    101                int iColumn;
    102                int numberColumns = rhs.quadraticObjective_->getNumCols();
    103                int numberBelow = 0;
    104                int numberAbove = 0;
    105                int numberDiagonal = 0;
    106                CoinZeroN(columnQuadraticLength2, numberExtendedColumns_);
    107                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    108                     for (CoinBigIndex j = columnQuadraticStart1[iColumn];
    109                               j < columnQuadraticStart1[iColumn] + columnQuadraticLength1[iColumn]; j++) {
    110                          int jColumn = columnQuadratic1[j];
    111                          if (jColumn > iColumn) {
    112                               numberBelow++;
    113                               columnQuadraticLength2[jColumn]++;
    114                               columnQuadraticLength2[iColumn]++;
    115                          } else if (jColumn == iColumn) {
    116                               numberDiagonal++;
    117                               columnQuadraticLength2[iColumn]++;
    118                          } else {
    119                               numberAbove++;
    120                          }
    121                     }
    122                }
    123                if (numberAbove > 0) {
    124                     if (numberAbove == numberBelow) {
    125                          // already done
    126                          quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
    127                          delete [] columnQuadraticStart2;
    128                          delete [] columnQuadraticLength2;
    129                     } else {
    130                          printf("number above = %d, number below = %d, error\n",
    131                                 numberAbove, numberBelow);
    132                          abort();
    133                     }
    134                } else {
    135                     int numberElements = numberDiagonal + 2 * numberBelow;
    136                     int * columnQuadratic2 = new int [numberElements];
    137                     double * quadraticElement2 = new double [numberElements];
    138                     columnQuadraticStart2[0] = 0;
    139                     numberElements = 0;
    140                     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    141                          int n = columnQuadraticLength2[iColumn];
    142                          columnQuadraticLength2[iColumn] = 0;
    143                          numberElements += n;
    144                          columnQuadraticStart2[iColumn+1] = numberElements;
    145                     }
    146                     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    147                          for (CoinBigIndex j = columnQuadraticStart1[iColumn];
    148                                    j < columnQuadraticStart1[iColumn] + columnQuadraticLength1[iColumn]; j++) {
    149                               int jColumn = columnQuadratic1[j];
    150                               if (jColumn > iColumn) {
    151                                    // put in two places
    152                                    CoinBigIndex put = columnQuadraticLength2[jColumn] + columnQuadraticStart2[jColumn];
    153                                    columnQuadraticLength2[jColumn]++;
    154                                    quadraticElement2[put] = quadraticElement1[j];
    155                                    columnQuadratic2[put] = iColumn;
    156                                    put = columnQuadraticLength2[iColumn] + columnQuadraticStart2[iColumn];
    157                                    columnQuadraticLength2[iColumn]++;
    158                                    quadraticElement2[put] = quadraticElement1[j];
    159                                    columnQuadratic2[put] = jColumn;
    160                               } else if (jColumn == iColumn) {
    161                                    CoinBigIndex put = columnQuadraticLength2[iColumn] + columnQuadraticStart2[iColumn];
    162                                    columnQuadraticLength2[iColumn]++;
    163                                    quadraticElement2[put] = quadraticElement1[j];
    164                                    columnQuadratic2[put] = iColumn;
    165                               } else {
    166                                    abort();
    167                               }
    168                          }
    169                     }
    170                     // Now create
    171                     quadraticObjective_ =
    172                          new CoinPackedMatrix (true,
    173                                                rhs.numberExtendedColumns_,
    174                                                rhs.numberExtendedColumns_,
    175                                                numberElements,
    176                                                quadraticElement2,
    177                                                columnQuadratic2,
    178                                                columnQuadraticStart2,
    179                                                columnQuadraticLength2, 0.0, 0.0);
    180                     delete [] columnQuadraticStart2;
    181                     delete [] columnQuadraticLength2;
    182                     delete [] columnQuadratic2;
    183                     delete [] quadraticElement2;
    184                }
     68ClpQuadraticObjective::ClpQuadraticObjective(const ClpQuadraticObjective &rhs,
     69  int type)
     70  : ClpObjective(rhs)
     71{
     72  numberColumns_ = rhs.numberColumns_;
     73  numberExtendedColumns_ = rhs.numberExtendedColumns_;
     74  fullMatrix_ = rhs.fullMatrix_;
     75  if (rhs.objective_) {
     76    objective_ = new double[numberExtendedColumns_];
     77    CoinMemcpyN(rhs.objective_, numberExtendedColumns_, objective_);
     78  } else {
     79    objective_ = NULL;
     80  }
     81  if (rhs.gradient_) {
     82    gradient_ = new double[numberExtendedColumns_];
     83    CoinMemcpyN(rhs.gradient_, numberExtendedColumns_, gradient_);
     84  } else {
     85    gradient_ = NULL;
     86  }
     87  if (rhs.quadraticObjective_) {
     88    // see what type of matrix wanted
     89    if (type == 0) {
     90      // just copy
     91      quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
     92    } else if (type == 1) {
     93      // expand to full symmetric
     94      fullMatrix_ = true;
     95      const int *columnQuadratic1 = rhs.quadraticObjective_->getIndices();
     96      const CoinBigIndex *columnQuadraticStart1 = rhs.quadraticObjective_->getVectorStarts();
     97      const int *columnQuadraticLength1 = rhs.quadraticObjective_->getVectorLengths();
     98      const double *quadraticElement1 = rhs.quadraticObjective_->getElements();
     99      CoinBigIndex *columnQuadraticStart2 = new CoinBigIndex[numberExtendedColumns_ + 1];
     100      int *columnQuadraticLength2 = new int[numberExtendedColumns_];
     101      int iColumn;
     102      int numberColumns = rhs.quadraticObjective_->getNumCols();
     103      int numberBelow = 0;
     104      int numberAbove = 0;
     105      int numberDiagonal = 0;
     106      CoinZeroN(columnQuadraticLength2, numberExtendedColumns_);
     107      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     108        for (CoinBigIndex j = columnQuadraticStart1[iColumn];
     109             j < columnQuadraticStart1[iColumn] + columnQuadraticLength1[iColumn]; j++) {
     110          int jColumn = columnQuadratic1[j];
     111          if (jColumn > iColumn) {
     112            numberBelow++;
     113            columnQuadraticLength2[jColumn]++;
     114            columnQuadraticLength2[iColumn]++;
     115          } else if (jColumn == iColumn) {
     116            numberDiagonal++;
     117            columnQuadraticLength2[iColumn]++;
    185118          } else {
    186                fullMatrix_ = false;
    187                abort(); // code when needed
    188           }
    189 
    190      } else {
    191           quadraticObjective_ = NULL;
    192      }
     119            numberAbove++;
     120          }
     121        }
     122      }
     123      if (numberAbove > 0) {
     124        if (numberAbove == numberBelow) {
     125          // already done
     126          quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
     127          delete[] columnQuadraticStart2;
     128          delete[] columnQuadraticLength2;
     129        } else {
     130          printf("number above = %d, number below = %d, error\n",
     131            numberAbove, numberBelow);
     132          abort();
     133        }
     134      } else {
     135        int numberElements = numberDiagonal + 2 * numberBelow;
     136        int *columnQuadratic2 = new int[numberElements];
     137        double *quadraticElement2 = new double[numberElements];
     138        columnQuadraticStart2[0] = 0;
     139        numberElements = 0;
     140        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     141          int n = columnQuadraticLength2[iColumn];
     142          columnQuadraticLength2[iColumn] = 0;
     143          numberElements += n;
     144          columnQuadraticStart2[iColumn + 1] = numberElements;
     145        }
     146        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     147          for (CoinBigIndex j = columnQuadraticStart1[iColumn];
     148               j < columnQuadraticStart1[iColumn] + columnQuadraticLength1[iColumn]; j++) {
     149            int jColumn = columnQuadratic1[j];
     150            if (jColumn > iColumn) {
     151              // put in two places
     152              CoinBigIndex put = columnQuadraticLength2[jColumn] + columnQuadraticStart2[jColumn];
     153              columnQuadraticLength2[jColumn]++;
     154              quadraticElement2[put] = quadraticElement1[j];
     155              columnQuadratic2[put] = iColumn;
     156              put = columnQuadraticLength2[iColumn] + columnQuadraticStart2[iColumn];
     157              columnQuadraticLength2[iColumn]++;
     158              quadraticElement2[put] = quadraticElement1[j];
     159              columnQuadratic2[put] = jColumn;
     160            } else if (jColumn == iColumn) {
     161              CoinBigIndex put = columnQuadraticLength2[iColumn] + columnQuadraticStart2[iColumn];
     162              columnQuadraticLength2[iColumn]++;
     163              quadraticElement2[put] = quadraticElement1[j];
     164              columnQuadratic2[put] = iColumn;
     165            } else {
     166              abort();
     167            }
     168          }
     169        }
     170        // Now create
     171        quadraticObjective_ = new CoinPackedMatrix(true,
     172          rhs.numberExtendedColumns_,
     173          rhs.numberExtendedColumns_,
     174          numberElements,
     175          quadraticElement2,
     176          columnQuadratic2,
     177          columnQuadraticStart2,
     178          columnQuadraticLength2, 0.0, 0.0);
     179        delete[] columnQuadraticStart2;
     180        delete[] columnQuadraticLength2;
     181        delete[] columnQuadratic2;
     182        delete[] quadraticElement2;
     183      }
     184    } else {
     185      fullMatrix_ = false;
     186      abort(); // code when needed
     187    }
     188
     189  } else {
     190    quadraticObjective_ = NULL;
     191  }
    193192}
    194193/* Subset constructor.  Duplicates are allowed
    195194   and order is as given.
    196195*/
    197 ClpQuadraticObjective::ClpQuadraticObjective (const ClpQuadraticObjective &rhs,
    198           int numberColumns,
    199           const int * whichColumn)
    200      : ClpObjective(rhs)
    201 {
    202      fullMatrix_ = rhs.fullMatrix_;
    203      objective_ = NULL;
    204      int extra = rhs.numberExtendedColumns_ - rhs.numberColumns_;
    205      numberColumns_ = 0;
    206      numberExtendedColumns_ = numberColumns + extra;
    207      if (numberColumns > 0) {
    208           // check valid lists
    209           int numberBad = 0;
    210           int i;
    211           for (i = 0; i < numberColumns; i++)
    212                if (whichColumn[i] < 0 || whichColumn[i] >= rhs.numberColumns_)
    213                     numberBad++;
    214           if (numberBad)
    215                throw CoinError("bad column list", "subset constructor",
    216                                "ClpQuadraticObjective");
    217           numberColumns_ = numberColumns;
    218           objective_ = new double[numberExtendedColumns_];
    219           for (i = 0; i < numberColumns_; i++)
    220                objective_[i] = rhs.objective_[whichColumn[i]];
    221           CoinMemcpyN(rhs.objective_ + rhs.numberColumns_,      (numberExtendedColumns_ - numberColumns_),
    222                       objective_ + numberColumns_);
    223           if (rhs.gradient_) {
    224                gradient_ = new double[numberExtendedColumns_];
    225                for (i = 0; i < numberColumns_; i++)
    226                     gradient_[i] = rhs.gradient_[whichColumn[i]];
    227                CoinMemcpyN(rhs.gradient_ + rhs.numberColumns_,  (numberExtendedColumns_ - numberColumns_),
    228                            gradient_ + numberColumns_);
    229           } else {
    230                gradient_ = NULL;
    231           }
    232      } else {
    233           gradient_ = NULL;
    234           objective_ = NULL;
    235      }
    236      if (rhs.quadraticObjective_) {
    237           quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_,
    238                     numberColumns, whichColumn,
    239                     numberColumns, whichColumn);
    240      } else {
    241           quadraticObjective_ = NULL;
    242      }
    243 }
    244 
     196ClpQuadraticObjective::ClpQuadraticObjective(const ClpQuadraticObjective &rhs,
     197  int numberColumns,
     198  const int *whichColumn)
     199  : ClpObjective(rhs)
     200{
     201  fullMatrix_ = rhs.fullMatrix_;
     202  objective_ = NULL;
     203  int extra = rhs.numberExtendedColumns_ - rhs.numberColumns_;
     204  numberColumns_ = 0;
     205  numberExtendedColumns_ = numberColumns + extra;
     206  if (numberColumns > 0) {
     207    // check valid lists
     208    int numberBad = 0;
     209    int i;
     210    for (i = 0; i < numberColumns; i++)
     211      if (whichColumn[i] < 0 || whichColumn[i] >= rhs.numberColumns_)
     212        numberBad++;
     213    if (numberBad)
     214      throw CoinError("bad column list", "subset constructor",
     215        "ClpQuadraticObjective");
     216    numberColumns_ = numberColumns;
     217    objective_ = new double[numberExtendedColumns_];
     218    for (i = 0; i < numberColumns_; i++)
     219      objective_[i] = rhs.objective_[whichColumn[i]];
     220    CoinMemcpyN(rhs.objective_ + rhs.numberColumns_, (numberExtendedColumns_ - numberColumns_),
     221      objective_ + numberColumns_);
     222    if (rhs.gradient_) {
     223      gradient_ = new double[numberExtendedColumns_];
     224      for (i = 0; i < numberColumns_; i++)
     225        gradient_[i] = rhs.gradient_[whichColumn[i]];
     226      CoinMemcpyN(rhs.gradient_ + rhs.numberColumns_, (numberExtendedColumns_ - numberColumns_),
     227        gradient_ + numberColumns_);
     228    } else {
     229      gradient_ = NULL;
     230    }
     231  } else {
     232    gradient_ = NULL;
     233    objective_ = NULL;
     234  }
     235  if (rhs.quadraticObjective_) {
     236    quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_,
     237      numberColumns, whichColumn,
     238      numberColumns, whichColumn);
     239  } else {
     240    quadraticObjective_ = NULL;
     241  }
     242}
    245243
    246244//-------------------------------------------------------------------
    247245// Destructor
    248246//-------------------------------------------------------------------
    249 ClpQuadraticObjective::~ClpQuadraticObjective ()
    250 {
    251      delete [] objective_;
    252      delete [] gradient_;
    253      delete quadraticObjective_;
     247ClpQuadraticObjective::~ClpQuadraticObjective()
     248{
     249  delete[] objective_;
     250  delete[] gradient_;
     251  delete quadraticObjective_;
    254252}
    255253
     
    258256//-------------------------------------------------------------------
    259257ClpQuadraticObjective &
    260 ClpQuadraticObjective::operator=(const ClpQuadraticObjective& rhs)
    261 {
    262      if (this != &rhs) {
    263           fullMatrix_ = rhs.fullMatrix_;
    264           delete quadraticObjective_;
    265           quadraticObjective_ = NULL;
    266           delete [] objective_;
    267           delete [] gradient_;
    268           ClpObjective::operator=(rhs);
    269           numberColumns_ = rhs.numberColumns_;
    270           numberExtendedColumns_ = rhs.numberExtendedColumns_;
    271           if (rhs.objective_) {
    272                objective_ = new double [numberExtendedColumns_];
    273                CoinMemcpyN(rhs.objective_, numberExtendedColumns_, objective_);
    274           } else {
    275                objective_ = NULL;
    276           }
    277           if (rhs.gradient_) {
    278                gradient_ = new double [numberExtendedColumns_];
    279                CoinMemcpyN(rhs.gradient_, numberExtendedColumns_, gradient_);
    280           } else {
    281                gradient_ = NULL;
    282           }
    283           if (rhs.quadraticObjective_) {
    284                quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
    285           } else {
    286                quadraticObjective_ = NULL;
    287           }
    288      }
    289      return *this;
     258ClpQuadraticObjective::operator=(const ClpQuadraticObjective &rhs)
     259{
     260  if (this != &rhs) {
     261    fullMatrix_ = rhs.fullMatrix_;
     262    delete quadraticObjective_;
     263    quadraticObjective_ = NULL;
     264    delete[] objective_;
     265    delete[] gradient_;
     266    ClpObjective::operator=(rhs);
     267    numberColumns_ = rhs.numberColumns_;
     268    numberExtendedColumns_ = rhs.numberExtendedColumns_;
     269    if (rhs.objective_) {
     270      objective_ = new double[numberExtendedColumns_];
     271      CoinMemcpyN(rhs.objective_, numberExtendedColumns_, objective_);
     272    } else {
     273      objective_ = NULL;
     274    }
     275    if (rhs.gradient_) {
     276      gradient_ = new double[numberExtendedColumns_];
     277      CoinMemcpyN(rhs.gradient_, numberExtendedColumns_, gradient_);
     278    } else {
     279      gradient_ = NULL;
     280    }
     281    if (rhs.quadraticObjective_) {
     282      quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
     283    } else {
     284      quadraticObjective_ = NULL;
     285    }
     286  }
     287  return *this;
    290288}
    291289
    292290// Returns gradient
    293291double *
    294 ClpQuadraticObjective::gradient(const ClpSimplex * model,
    295                                 const double * solution, double & offset, bool refresh,
    296                                 int includeLinear)
    297 {
    298      offset = 0.0;
    299      bool scaling = false;
    300      if (model && (model->rowScale() ||
    301                    model->objectiveScale() != 1.0 || model->optimizationDirection() != 1.0))
    302           scaling = true;
    303      const double * cost = NULL;
    304      if (model)
    305           cost = model->costRegion();
    306      if (!cost) {
    307           // not in solve
    308           cost = objective_;
    309           scaling = false;
    310      }
    311      if (!scaling) {
    312           if (!quadraticObjective_ || !solution || !activated_) {
    313                return objective_;
     292ClpQuadraticObjective::gradient(const ClpSimplex *model,
     293  const double *solution, double &offset, bool refresh,
     294  int includeLinear)
     295{
     296  offset = 0.0;
     297  bool scaling = false;
     298  if (model && (model->rowScale() || model->objectiveScale() != 1.0 || model->optimizationDirection() != 1.0))
     299    scaling = true;
     300  const double *cost = NULL;
     301  if (model)
     302    cost = model->costRegion();
     303  if (!cost) {
     304    // not in solve
     305    cost = objective_;
     306    scaling = false;
     307  }
     308  if (!scaling) {
     309    if (!quadraticObjective_ || !solution || !activated_) {
     310      return objective_;
     311    } else {
     312      if (refresh || !gradient_) {
     313        if (!gradient_)
     314          gradient_ = new double[numberExtendedColumns_];
     315        const int *columnQuadratic = quadraticObjective_->getIndices();
     316        const CoinBigIndex *columnQuadraticStart = quadraticObjective_->getVectorStarts();
     317        const int *columnQuadraticLength = quadraticObjective_->getVectorLengths();
     318        const double *quadraticElement = quadraticObjective_->getElements();
     319        offset = 0.0;
     320        // use current linear cost region
     321        if (includeLinear == 1)
     322          CoinMemcpyN(cost, numberExtendedColumns_, gradient_);
     323        else if (includeLinear == 2)
     324          CoinMemcpyN(objective_, numberExtendedColumns_, gradient_);
     325        else
     326          memset(gradient_, 0, numberExtendedColumns_ * sizeof(double));
     327        if (activated_) {
     328          if (!fullMatrix_) {
     329            int iColumn;
     330            for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     331              double valueI = solution[iColumn];
     332              CoinBigIndex j;
     333              for (j = columnQuadraticStart[iColumn];
     334                   j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     335                int jColumn = columnQuadratic[j];
     336                double valueJ = solution[jColumn];
     337                double elementValue = quadraticElement[j];
     338                if (iColumn != jColumn) {
     339                  offset += valueI * valueJ * elementValue;
     340                  //if (fabs(valueI*valueJ*elementValue)>1.0e-12)
     341                  //printf("%d %d %g %g %g -> %g\n",
     342                  //       iColumn,jColumn,valueI,valueJ,elementValue,
     343                  //       valueI*valueJ*elementValue);
     344                  double gradientI = valueJ * elementValue;
     345                  double gradientJ = valueI * elementValue;
     346                  gradient_[iColumn] += gradientI;
     347                  gradient_[jColumn] += gradientJ;
     348                } else {
     349                  offset += 0.5 * valueI * valueI * elementValue;
     350                  //if (fabs(valueI*valueI*elementValue)>1.0e-12)
     351                  //printf("XX %d %g %g -> %g\n",
     352                  //       iColumn,valueI,elementValue,
     353                  //       0.5*valueI*valueI*elementValue);
     354                  double gradientI = valueI * elementValue;
     355                  gradient_[iColumn] += gradientI;
     356                }
     357              }
     358            }
    314359          } else {
    315                if (refresh || !gradient_) {
    316                     if (!gradient_)
    317                          gradient_ = new double[numberExtendedColumns_];
    318                     const int * columnQuadratic = quadraticObjective_->getIndices();
    319                     const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
    320                     const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
    321                     const double * quadraticElement = quadraticObjective_->getElements();
    322                     offset = 0.0;
    323                     // use current linear cost region
    324                     if (includeLinear == 1)
    325                          CoinMemcpyN(cost, numberExtendedColumns_, gradient_);
    326                     else if (includeLinear == 2)
    327                          CoinMemcpyN(objective_, numberExtendedColumns_, gradient_);
    328                     else
    329                          memset(gradient_, 0, numberExtendedColumns_ * sizeof(double));
    330                     if (activated_) {
    331                          if (!fullMatrix_) {
    332                               int iColumn;
    333                               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    334                                    double valueI = solution[iColumn];
    335                                    CoinBigIndex j;
    336                                    for (j = columnQuadraticStart[iColumn];
    337                                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    338                                         int jColumn = columnQuadratic[j];
    339                                         double valueJ = solution[jColumn];
    340                                         double elementValue = quadraticElement[j];
    341                                         if (iColumn != jColumn) {
    342                                              offset += valueI * valueJ * elementValue;
    343                                              //if (fabs(valueI*valueJ*elementValue)>1.0e-12)
    344                                              //printf("%d %d %g %g %g -> %g\n",
    345                                              //       iColumn,jColumn,valueI,valueJ,elementValue,
    346                                              //       valueI*valueJ*elementValue);
    347                                              double gradientI = valueJ * elementValue;
    348                                              double gradientJ = valueI * elementValue;
    349                                              gradient_[iColumn] += gradientI;
    350                                              gradient_[jColumn] += gradientJ;
    351                                         } else {
    352                                              offset += 0.5 * valueI * valueI * elementValue;
    353                                              //if (fabs(valueI*valueI*elementValue)>1.0e-12)
    354                                              //printf("XX %d %g %g -> %g\n",
    355                                              //       iColumn,valueI,elementValue,
    356                                              //       0.5*valueI*valueI*elementValue);
    357                                              double gradientI = valueI * elementValue;
    358                                              gradient_[iColumn] += gradientI;
    359                                         }
    360                                    }
    361                               }
    362                          } else {
    363                               // full matrix
    364                               int iColumn;
    365                               offset *= 2.0;
    366                               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    367                                    CoinBigIndex j;
    368                                    double value = 0.0;
    369                                    double current = gradient_[iColumn];
    370                                    for (j = columnQuadraticStart[iColumn];
    371                                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    372                                         int jColumn = columnQuadratic[j];
    373                                         double valueJ = solution[jColumn] * quadraticElement[j];
    374                                         value += valueJ;
    375                                    }
    376                                    offset += value * solution[iColumn];
    377                                    gradient_[iColumn] = current + value;
    378                               }
    379                               offset *= 0.5;
    380                          }
    381                     }
    382                }
    383                if (model)
    384                     offset *= model->optimizationDirection() * model->objectiveScale();
    385                return gradient_;
    386           }
    387      } else {
    388           // do scaling
    389           assert (solution);
    390           // for now only if half
    391           assert (!fullMatrix_);
    392           if (refresh || !gradient_) {
    393                if (!gradient_)
    394                     gradient_ = new double[numberExtendedColumns_];
    395                double direction = model->optimizationDirection() * model->objectiveScale();
    396                // direction is actually scale out not scale in
    397                //if (direction)
    398                //direction = 1.0/direction;
    399                const int * columnQuadratic = quadraticObjective_->getIndices();
    400                const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
    401                const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
    402                const double * quadraticElement = quadraticObjective_->getElements();
    403                int iColumn;
    404                const double * columnScale = model->columnScale();
    405                // use current linear cost region (already scaled)
    406                if (includeLinear == 1) {
    407                     CoinMemcpyN(model->costRegion(), numberExtendedColumns_, gradient_);
    408                }        else if (includeLinear == 2) {
    409                     memset(gradient_ + numberColumns_, 0, (numberExtendedColumns_ - numberColumns_)*sizeof(double));
    410                     if (!columnScale) {
    411                          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    412                               gradient_[iColumn] = objective_[iColumn] * direction;
    413                          }
    414                     } else {
    415                          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    416                               gradient_[iColumn] = objective_[iColumn] * direction * columnScale[iColumn];
    417                          }
    418                     }
    419                } else {
    420                     memset(gradient_, 0, numberExtendedColumns_ * sizeof(double));
    421                }
    422                if (!columnScale) {
    423                     if (activated_) {
    424                          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    425                               double valueI = solution[iColumn];
    426                               CoinBigIndex j;
    427                               for (j = columnQuadraticStart[iColumn];
    428                                         j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    429                                    int jColumn = columnQuadratic[j];
    430                                    double valueJ = solution[jColumn];
    431                                    double elementValue = quadraticElement[j];
    432                                    elementValue *= direction;
    433                                    if (iColumn != jColumn) {
    434                                         offset += valueI * valueJ * elementValue;
    435                                         double gradientI = valueJ * elementValue;
    436                                         double gradientJ = valueI * elementValue;
    437                                         gradient_[iColumn] += gradientI;
    438                                         gradient_[jColumn] += gradientJ;
    439                                    } else {
    440                                         offset += 0.5 * valueI * valueI * elementValue;
    441                                         double gradientI = valueI * elementValue;
    442                                         gradient_[iColumn] += gradientI;
    443                                    }
    444                               }
    445                          }
    446                     }
    447                } else {
    448                     // scaling
    449                     if (activated_) {
    450                          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    451                               double valueI = solution[iColumn];
    452                               double scaleI = columnScale[iColumn] * direction;
    453                               CoinBigIndex j;
    454                               for (j = columnQuadraticStart[iColumn];
    455                                         j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    456                                    int jColumn = columnQuadratic[j];
    457                                    double valueJ = solution[jColumn];
    458                                    double elementValue = quadraticElement[j];
    459                                    double scaleJ = columnScale[jColumn];
    460                                    elementValue *= scaleI * scaleJ;
    461                                    if (iColumn != jColumn) {
    462                                         offset += valueI * valueJ * elementValue;
    463                                         double gradientI = valueJ * elementValue;
    464                                         double gradientJ = valueI * elementValue;
    465                                         gradient_[iColumn] += gradientI;
    466                                         gradient_[jColumn] += gradientJ;
    467                                    } else {
    468                                         offset += 0.5 * valueI * valueI * elementValue;
    469                                         double gradientI = valueI * elementValue;
    470                                         gradient_[iColumn] += gradientI;
    471                                    }
    472                               }
    473                          }
    474                     }
    475                }
    476           }
    477           if (model)
    478                offset *= model->optimizationDirection();
    479           return gradient_;
    480      }
     360            // full matrix
     361            int iColumn;
     362            offset *= 2.0;
     363            for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     364              CoinBigIndex j;
     365              double value = 0.0;
     366              double current = gradient_[iColumn];
     367              for (j = columnQuadraticStart[iColumn];
     368                   j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     369                int jColumn = columnQuadratic[j];
     370                double valueJ = solution[jColumn] * quadraticElement[j];
     371                value += valueJ;
     372              }
     373              offset += value * solution[iColumn];
     374              gradient_[iColumn] = current + value;
     375            }
     376            offset *= 0.5;
     377          }
     378        }
     379      }
     380      if (model)
     381        offset *= model->optimizationDirection() * model->objectiveScale();
     382      return gradient_;
     383    }
     384  } else {
     385    // do scaling
     386    assert(solution);
     387    // for now only if half
     388    assert(!fullMatrix_);
     389    if (refresh || !gradient_) {
     390      if (!gradient_)
     391        gradient_ = new double[numberExtendedColumns_];
     392      double direction = model->optimizationDirection() * model->objectiveScale();
     393      // direction is actually scale out not scale in
     394      //if (direction)
     395      //direction = 1.0/direction;
     396      const int *columnQuadratic = quadraticObjective_->getIndices();
     397      const CoinBigIndex *columnQuadraticStart = quadraticObjective_->getVectorStarts();
     398      const int *columnQuadraticLength = quadraticObjective_->getVectorLengths();
     399      const double *quadraticElement = quadraticObjective_->getElements();
     400      int iColumn;
     401      const double *columnScale = model->columnScale();
     402      // use current linear cost region (already scaled)
     403      if (includeLinear == 1) {
     404        CoinMemcpyN(model->costRegion(), numberExtendedColumns_, gradient_);
     405      } else if (includeLinear == 2) {
     406        memset(gradient_ + numberColumns_, 0, (numberExtendedColumns_ - numberColumns_) * sizeof(double));
     407        if (!columnScale) {
     408          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     409            gradient_[iColumn] = objective_[iColumn] * direction;
     410          }
     411        } else {
     412          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     413            gradient_[iColumn] = objective_[iColumn] * direction * columnScale[iColumn];
     414          }
     415        }
     416      } else {
     417        memset(gradient_, 0, numberExtendedColumns_ * sizeof(double));
     418      }
     419      if (!columnScale) {
     420        if (activated_) {
     421          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     422            double valueI = solution[iColumn];
     423            CoinBigIndex j;
     424            for (j = columnQuadraticStart[iColumn];
     425                 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     426              int jColumn = columnQuadratic[j];
     427              double valueJ = solution[jColumn];
     428              double elementValue = quadraticElement[j];
     429              elementValue *= direction;
     430              if (iColumn != jColumn) {
     431                offset += valueI * valueJ * elementValue;
     432                double gradientI = valueJ * elementValue;
     433                double gradientJ = valueI * elementValue;
     434                gradient_[iColumn] += gradientI;
     435                gradient_[jColumn] += gradientJ;
     436              } else {
     437                offset += 0.5 * valueI * valueI * elementValue;
     438                double gradientI = valueI * elementValue;
     439                gradient_[iColumn] += gradientI;
     440              }
     441            }
     442          }
     443        }
     444      } else {
     445        // scaling
     446        if (activated_) {
     447          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     448            double valueI = solution[iColumn];
     449            double scaleI = columnScale[iColumn] * direction;
     450            CoinBigIndex j;
     451            for (j = columnQuadraticStart[iColumn];
     452                 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     453              int jColumn = columnQuadratic[j];
     454              double valueJ = solution[jColumn];
     455              double elementValue = quadraticElement[j];
     456              double scaleJ = columnScale[jColumn];
     457              elementValue *= scaleI * scaleJ;
     458              if (iColumn != jColumn) {
     459                offset += valueI * valueJ * elementValue;
     460                double gradientI = valueJ * elementValue;
     461                double gradientJ = valueI * elementValue;
     462                gradient_[iColumn] += gradientI;
     463                gradient_[jColumn] += gradientJ;
     464              } else {
     465                offset += 0.5 * valueI * valueI * elementValue;
     466                double gradientI = valueI * elementValue;
     467                gradient_[iColumn] += gradientI;
     468              }
     469            }
     470          }
     471        }
     472      }
     473    }
     474    if (model)
     475      offset *= model->optimizationDirection();
     476    return gradient_;
     477  }
    481478}
    482479
     
    484481// Clone
    485482//-------------------------------------------------------------------
    486 ClpObjective * ClpQuadraticObjective::clone() const
    487 {
    488      return new ClpQuadraticObjective(*this);
     483ClpObjective *ClpQuadraticObjective::clone() const
     484{
     485  return new ClpQuadraticObjective(*this);
    489486}
    490487/* Subset clone.  Duplicates are allowed
     
    492489*/
    493490ClpObjective *
    494 ClpQuadraticObjective::subsetClone (int numberColumns,
    495                                     const int * whichColumns) const
    496 {
    497      return new ClpQuadraticObjective(*this, numberColumns, whichColumns);
     491ClpQuadraticObjective::subsetClone(int numberColumns,
     492  const int *whichColumns) const
     493{
     494  return new ClpQuadraticObjective(*this, numberColumns, whichColumns);
    498495}
    499496// Resize objective
    500 void
    501 ClpQuadraticObjective::resize(int newNumberColumns)
    502 {
    503      if (numberColumns_ != newNumberColumns) {
    504           int newExtended = newNumberColumns + (numberExtendedColumns_ - numberColumns_);
    505           int i;
    506           double * newArray = new double[newExtended];
    507           if (objective_)
    508                CoinMemcpyN(objective_,  CoinMin(newExtended, numberExtendedColumns_), newArray);
    509           delete [] objective_;
    510           objective_ = newArray;
    511           for (i = numberColumns_; i < newNumberColumns; i++)
    512                objective_[i] = 0.0;
    513           if (gradient_) {
    514                newArray = new double[newExtended];
    515                if (gradient_)
    516                     CoinMemcpyN(gradient_,      CoinMin(newExtended, numberExtendedColumns_), newArray);
    517                delete [] gradient_;
    518                gradient_ = newArray;
    519                for (i = numberColumns_; i < newNumberColumns; i++)
    520                     gradient_[i] = 0.0;
    521           }
    522           if (quadraticObjective_) {
    523                if (newNumberColumns < numberColumns_) {
    524                     int * which = new int[numberColumns_-newNumberColumns];
    525                     int i;
    526                     for (i = newNumberColumns; i < numberColumns_; i++)
    527                          which[i-newNumberColumns] = i;
    528                     quadraticObjective_->deleteRows(numberColumns_ - newNumberColumns, which);
    529                     quadraticObjective_->deleteCols(numberColumns_ - newNumberColumns, which);
    530                     delete [] which;
    531                } else {
    532                     quadraticObjective_->setDimensions(newNumberColumns, newNumberColumns);
    533                }
    534           }
    535           numberColumns_ = newNumberColumns;
    536           numberExtendedColumns_ = newExtended;
    537      }
    538 
     497void ClpQuadraticObjective::resize(int newNumberColumns)
     498{
     499  if (numberColumns_ != newNumberColumns) {
     500    int newExtended = newNumberColumns + (numberExtendedColumns_ - numberColumns_);
     501    int i;
     502    double *newArray = new double[newExtended];
     503    if (objective_)
     504      CoinMemcpyN(objective_, CoinMin(newExtended, numberExtendedColumns_), newArray);
     505    delete[] objective_;
     506    objective_ = newArray;
     507    for (i = numberColumns_; i < newNumberColumns; i++)
     508      objective_[i] = 0.0;
     509    if (gradient_) {
     510      newArray = new double[newExtended];
     511      if (gradient_)
     512        CoinMemcpyN(gradient_, CoinMin(newExtended, numberExtendedColumns_), newArray);
     513      delete[] gradient_;
     514      gradient_ = newArray;
     515      for (i = numberColumns_; i < newNumberColumns; i++)
     516        gradient_[i] = 0.0;
     517    }
     518    if (quadraticObjective_) {
     519      if (newNumberColumns < numberColumns_) {
     520        int *which = new int[numberColumns_ - newNumberColumns];
     521        int i;
     522        for (i = newNumberColumns; i < numberColumns_; i++)
     523          which[i - newNumberColumns] = i;
     524        quadraticObjective_->deleteRows(numberColumns_ - newNumberColumns, which);
     525        quadraticObjective_->deleteCols(numberColumns_ - newNumberColumns, which);
     526        delete[] which;
     527      } else {
     528        quadraticObjective_->setDimensions(newNumberColumns, newNumberColumns);
     529      }
     530    }
     531    numberColumns_ = newNumberColumns;
     532    numberExtendedColumns_ = newExtended;
     533  }
    539534}
    540535// Delete columns in  objective
    541 void
    542 ClpQuadraticObjective::deleteSome(int numberToDelete, const int * which)
    543 {
    544      int newNumberColumns = numberColumns_ - numberToDelete;
    545      int newExtended = numberExtendedColumns_ - numberToDelete;
    546      if (objective_) {
    547           int i ;
    548           char * deleted = new char[numberColumns_];
    549           int numberDeleted = 0;
    550           memset(deleted, 0, numberColumns_ * sizeof(char));
    551           for (i = 0; i < numberToDelete; i++) {
    552                int j = which[i];
    553                if (j >= 0 && j < numberColumns_ && !deleted[j]) {
    554                     numberDeleted++;
    555                     deleted[j] = 1;
    556                }
    557           }
    558           newNumberColumns = numberColumns_ - numberDeleted;
    559           newExtended = numberExtendedColumns_ - numberDeleted;
    560           double * newArray = new double[newExtended];
    561           int put = 0;
    562           for (i = 0; i < numberColumns_; i++) {
    563                if (!deleted[i]) {
    564                     newArray[put++] = objective_[i];
    565                }
    566           }
    567           delete [] objective_;
    568           objective_ = newArray;
    569           delete [] deleted;
    570           CoinMemcpyN(objective_ + numberColumns_,      (numberExtendedColumns_ - numberColumns_),
    571                       objective_ + newNumberColumns);
    572      }
    573      if (gradient_) {
    574           int i ;
    575           char * deleted = new char[numberColumns_];
    576           int numberDeleted = 0;
    577           memset(deleted, 0, numberColumns_ * sizeof(char));
    578           for (i = 0; i < numberToDelete; i++) {
    579                int j = which[i];
    580                if (j >= 0 && j < numberColumns_ && !deleted[j]) {
    581                     numberDeleted++;
    582                     deleted[j] = 1;
    583                }
    584           }
    585           newNumberColumns = numberColumns_ - numberDeleted;
    586           newExtended = numberExtendedColumns_ - numberDeleted;
    587           double * newArray = new double[newExtended];
    588           int put = 0;
    589           for (i = 0; i < numberColumns_; i++) {
    590                if (!deleted[i]) {
    591                     newArray[put++] = gradient_[i];
    592                }
    593           }
    594           delete [] gradient_;
    595           gradient_ = newArray;
    596           delete [] deleted;
    597           CoinMemcpyN(gradient_ + numberColumns_,       (numberExtendedColumns_ - numberColumns_),
    598                       gradient_ + newNumberColumns);
    599      }
    600      numberColumns_ = newNumberColumns;
    601      numberExtendedColumns_ = newExtended;
    602      if (quadraticObjective_) {
    603           quadraticObjective_->deleteCols(numberToDelete, which);
    604           quadraticObjective_->deleteRows(numberToDelete, which);
    605      }
     536void ClpQuadraticObjective::deleteSome(int numberToDelete, const int *which)
     537{
     538  int newNumberColumns = numberColumns_ - numberToDelete;
     539  int newExtended = numberExtendedColumns_ - numberToDelete;
     540  if (objective_) {
     541    int i;
     542    char *deleted = new char[numberColumns_];
     543    int numberDeleted = 0;
     544    memset(deleted, 0, numberColumns_ * sizeof(char));
     545    for (i = 0; i < numberToDelete; i++) {
     546      int j = which[i];
     547      if (j >= 0 && j < numberColumns_ && !deleted[j]) {
     548        numberDeleted++;
     549        deleted[j] = 1;
     550      }
     551    }
     552    newNumberColumns = numberColumns_ - numberDeleted;
     553    newExtended = numberExtendedColumns_ - numberDeleted;
     554    double *newArray = new double[newExtended];
     555    int put = 0;
     556    for (i = 0; i < numberColumns_; i++) {
     557      if (!deleted[i]) {
     558        newArray[put++] = objective_[i];
     559      }
     560    }
     561    delete[] objective_;
     562    objective_ = newArray;
     563    delete[] deleted;
     564    CoinMemcpyN(objective_ + numberColumns_, (numberExtendedColumns_ - numberColumns_),
     565      objective_ + newNumberColumns);
     566  }
     567  if (gradient_) {
     568    int i;
     569    char *deleted = new char[numberColumns_];
     570    int numberDeleted = 0;
     571    memset(deleted, 0, numberColumns_ * sizeof(char));
     572    for (i = 0; i < numberToDelete; i++) {
     573      int j = which[i];
     574      if (j >= 0 && j < numberColumns_ && !deleted[j]) {
     575        numberDeleted++;
     576        deleted[j] = 1;
     577      }
     578    }
     579    newNumberColumns = numberColumns_ - numberDeleted;
     580    newExtended = numberExtendedColumns_ - numberDeleted;
     581    double *newArray = new double[newExtended];
     582    int put = 0;
     583    for (i = 0; i < numberColumns_; i++) {
     584      if (!deleted[i]) {
     585        newArray[put++] = gradient_[i];
     586      }
     587    }
     588    delete[] gradient_;
     589    gradient_ = newArray;
     590    delete[] deleted;
     591    CoinMemcpyN(gradient_ + numberColumns_, (numberExtendedColumns_ - numberColumns_),
     592      gradient_ + newNumberColumns);
     593  }
     594  numberColumns_ = newNumberColumns;
     595  numberExtendedColumns_ = newExtended;
     596  if (quadraticObjective_) {
     597    quadraticObjective_->deleteCols(numberToDelete, which);
     598    quadraticObjective_->deleteRows(numberToDelete, which);
     599  }
    606600}
    607601
    608602// Load up quadratic objective
    609 void
    610 ClpQuadraticObjective::loadQuadraticObjective(const int numberColumns, const CoinBigIndex * start,
    611           const int * column, const double * element, int numberExtended)
    612 {
    613      fullMatrix_ = false;
    614      delete quadraticObjective_;
    615      quadraticObjective_ = new CoinPackedMatrix(true, numberColumns, numberColumns,
    616                start[numberColumns], element, column, start, NULL);
    617      numberColumns_ = numberColumns;
    618      if (numberExtended > numberExtendedColumns_) {
    619           if (objective_) {
    620                // make correct size
    621                double * newArray = new double[numberExtended];
    622                CoinMemcpyN(objective_, numberColumns_, newArray);
    623                delete [] objective_;
    624                objective_ = newArray;
    625                memset(objective_ + numberColumns_, 0, (numberExtended - numberColumns_)*sizeof(double));
    626           }
    627           if (gradient_) {
    628                // make correct size
    629                double * newArray = new double[numberExtended];
    630                CoinMemcpyN(gradient_, numberColumns_, newArray);
    631                delete [] gradient_;
    632                gradient_ = newArray;
    633                memset(gradient_ + numberColumns_, 0, (numberExtended - numberColumns_)*sizeof(double));
    634           }
    635           numberExtendedColumns_ = numberExtended;
    636      } else {
    637           numberExtendedColumns_ = numberColumns_;
    638      }
    639 }
    640 void
    641 ClpQuadraticObjective::loadQuadraticObjective (  const CoinPackedMatrix& matrix)
    642 {
    643      delete quadraticObjective_;
    644      quadraticObjective_ = new CoinPackedMatrix(matrix);
     603void ClpQuadraticObjective::loadQuadraticObjective(const int numberColumns, const CoinBigIndex *start,
     604  const int *column, const double *element, int numberExtended)
     605{
     606  fullMatrix_ = false;
     607  delete quadraticObjective_;
     608  quadraticObjective_ = new CoinPackedMatrix(true, numberColumns, numberColumns,
     609    start[numberColumns], element, column, start, NULL);
     610  numberColumns_ = numberColumns;
     611  if (numberExtended > numberExtendedColumns_) {
     612    if (objective_) {
     613      // make correct size
     614      double *newArray = new double[numberExtended];
     615      CoinMemcpyN(objective_, numberColumns_, newArray);
     616      delete[] objective_;
     617      objective_ = newArray;
     618      memset(objective_ + numberColumns_, 0, (numberExtended - numberColumns_) * sizeof(double));
     619    }
     620    if (gradient_) {
     621      // make correct size
     622      double *newArray = new double[numberExtended];
     623      CoinMemcpyN(gradient_, numberColumns_, newArray);
     624      delete[] gradient_;
     625      gradient_ = newArray;
     626      memset(gradient_ + numberColumns_, 0, (numberExtended - numberColumns_) * sizeof(double));
     627    }
     628    numberExtendedColumns_ = numberExtended;
     629  } else {
     630    numberExtendedColumns_ = numberColumns_;
     631  }
     632}
     633void ClpQuadraticObjective::loadQuadraticObjective(const CoinPackedMatrix &matrix)
     634{
     635  delete quadraticObjective_;
     636  quadraticObjective_ = new CoinPackedMatrix(matrix);
    645637}
    646638// Get rid of quadratic objective
    647 void
    648 ClpQuadraticObjective::deleteQuadraticObjective()
    649 {
    650      delete quadraticObjective_;
    651      quadraticObjective_ = NULL;
     639void ClpQuadraticObjective::deleteQuadraticObjective()
     640{
     641  delete quadraticObjective_;
     642  quadraticObjective_ = NULL;
    652643}
    653644/* Returns reduced gradient.Returns an offset (to be added to current one).
    654645 */
    655646double
    656 ClpQuadraticObjective::reducedGradient(ClpSimplex * model, double * region,
    657                                        bool useFeasibleCosts)
    658 {
    659      int numberRows = model->numberRows();
    660      int numberColumns = model->numberColumns();
    661 
    662      //work space
    663      CoinIndexedVector  * workSpace = model->rowArray(0);
    664 
    665      CoinIndexedVector arrayVector;
    666      arrayVector.reserve(numberRows + 1);
    667 
    668      int iRow;
     647ClpQuadraticObjective::reducedGradient(ClpSimplex *model, double *region,
     648  bool useFeasibleCosts)
     649{
     650  int numberRows = model->numberRows();
     651  int numberColumns = model->numberColumns();
     652
     653  //work space
     654  CoinIndexedVector *workSpace = model->rowArray(0);
     655
     656  CoinIndexedVector arrayVector;
     657  arrayVector.reserve(numberRows + 1);
     658
     659  int iRow;
    669660#ifdef CLP_DEBUG
    670      workSpace->checkClear();
     661  workSpace->checkClear();
    671662#endif
    672      double * array = arrayVector.denseVector();
    673      int * index = arrayVector.getIndices();
    674      int number = 0;
    675      const double * costNow = gradient(model, model->solutionRegion(), offset_,
    676                                        true, useFeasibleCosts ? 2 : 1);
    677      double * cost = model->costRegion();
    678      const int * pivotVariable = model->pivotVariable();
    679      for (iRow = 0; iRow < numberRows; iRow++) {
    680           int iPivot = pivotVariable[iRow];
    681           double value;
    682           if (iPivot < numberColumns)
    683                value = costNow[iPivot];
    684           else if (!useFeasibleCosts)
    685                value = cost[iPivot];
    686           else
    687                value = 0.0;
    688           if (value) {
    689                array[iRow] = value;
    690                index[number++] = iRow;
    691           }
    692      }
    693      arrayVector.setNumElements(number);
    694 
    695      // Btran basic costs
    696      model->factorization()->updateColumnTranspose(workSpace, &arrayVector);
    697      double * work = workSpace->denseVector();
    698      ClpFillN(work, numberRows, 0.0);
    699      // now look at dual solution
    700      double * rowReducedCost = region + numberColumns;
    701      double * dual = rowReducedCost;
    702      const double * rowCost = cost + numberColumns;
    703      for (iRow = 0; iRow < numberRows; iRow++) {
    704           dual[iRow] = array[iRow];
    705      }
    706      double * dj = region;
    707      ClpDisjointCopyN(costNow, numberColumns, dj);
    708 
    709      model->transposeTimes(-1.0, dual, dj);
    710      for (iRow = 0; iRow < numberRows; iRow++) {
    711           // slack
    712           double value = dual[iRow];
    713           value += rowCost[iRow];
    714           rowReducedCost[iRow] = value;
    715      }
    716      return offset_;
     663  double *array = arrayVector.denseVector();
     664  int *index = arrayVector.getIndices();
     665  int number = 0;
     666  const double *costNow = gradient(model, model->solutionRegion(), offset_,
     667    true, useFeasibleCosts ? 2 : 1);
     668  double *cost = model->costRegion();
     669  const int *pivotVariable = model->pivotVariable();
     670  for (iRow = 0; iRow < numberRows; iRow++) {
     671    int iPivot = pivotVariable[iRow];
     672    double value;
     673    if (iPivot < numberColumns)
     674      value = costNow[iPivot];
     675    else if (!useFeasibleCosts)
     676      value = cost[iPivot];
     677    else
     678      value = 0.0;
     679    if (value) {
     680      array[iRow] = value;
     681      index[number++] = iRow;
     682    }
     683  }
     684  arrayVector.setNumElements(number);
     685
     686  // Btran basic costs
     687  model->factorization()->updateColumnTranspose(workSpace, &arrayVector);
     688  double *work = workSpace->denseVector();
     689  ClpFillN(work, numberRows, 0.0);
     690  // now look at dual solution
     691  double *rowReducedCost = region + numberColumns;
     692  double *dual = rowReducedCost;
     693  const double *rowCost = cost + numberColumns;
     694  for (iRow = 0; iRow < numberRows; iRow++) {
     695    dual[iRow] = array[iRow];
     696  }
     697  double *dj = region;
     698  ClpDisjointCopyN(costNow, numberColumns, dj);
     699
     700  model->transposeTimes(-1.0, dual, dj);
     701  for (iRow = 0; iRow < numberRows; iRow++) {
     702    // slack
     703    double value = dual[iRow];
     704    value += rowCost[iRow];
     705    rowReducedCost[iRow] = value;
     706  }
     707  return offset_;
    717708}
    718709/* Returns step length which gives minimum of objective for
     
    722713*/
    723714double
    724 ClpQuadraticObjective::stepLength(ClpSimplex * model,
    725                                   const double * solution,
    726                                   const double * change,
    727                                   double maximumTheta,
    728                                   double & currentObj,
    729                                   double & predictedObj,
    730                                   double & thetaObj)
    731 {
    732      const double * cost = model->costRegion();
    733      bool inSolve = true;
    734      if (!cost) {
    735           // not in solve
    736           cost = objective_;
    737           inSolve = false;
    738      }
    739      double delta = 0.0;
    740      double linearCost = 0.0;
    741      int numberRows = model->numberRows();
    742      int numberColumns = model->numberColumns();
    743      int numberTotal = numberColumns;
    744      if (inSolve)
    745           numberTotal += numberRows;
    746      currentObj = 0.0;
    747      thetaObj = 0.0;
    748      for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
    749           delta += cost[iColumn] * change[iColumn];
    750           linearCost += cost[iColumn] * solution[iColumn];
    751      }
    752      if (!activated_ || !quadraticObjective_) {
    753           currentObj = linearCost;
    754           thetaObj = currentObj + delta * maximumTheta;
    755           if (delta < 0.0) {
    756                return maximumTheta;
     715ClpQuadraticObjective::stepLength(ClpSimplex *model,
     716  const double *solution,
     717  const double *change,
     718  double maximumTheta,
     719  double &currentObj,
     720  double &predictedObj,
     721  double &thetaObj)
     722{
     723  const double *cost = model->costRegion();
     724  bool inSolve = true;
     725  if (!cost) {
     726    // not in solve
     727    cost = objective_;
     728    inSolve = false;
     729  }
     730  double delta = 0.0;
     731  double linearCost = 0.0;
     732  int numberRows = model->numberRows();
     733  int numberColumns = model->numberColumns();
     734  int numberTotal = numberColumns;
     735  if (inSolve)
     736    numberTotal += numberRows;
     737  currentObj = 0.0;
     738  thetaObj = 0.0;
     739  for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
     740    delta += cost[iColumn] * change[iColumn];
     741    linearCost += cost[iColumn] * solution[iColumn];
     742  }
     743  if (!activated_ || !quadraticObjective_) {
     744    currentObj = linearCost;
     745    thetaObj = currentObj + delta * maximumTheta;
     746    if (delta < 0.0) {
     747      return maximumTheta;
     748    } else {
     749      COIN_DETAIL_PRINT(printf("odd linear direction %g\n", delta));
     750      return 0.0;
     751    }
     752  }
     753  assert(model);
     754  bool scaling = false;
     755  if ((model->rowScale() || model->objectiveScale() != 1.0 || model->optimizationDirection() != 1.0) && inSolve)
     756    scaling = true;
     757  const int *columnQuadratic = quadraticObjective_->getIndices();
     758  const CoinBigIndex *columnQuadraticStart = quadraticObjective_->getVectorStarts();
     759  const int *columnQuadraticLength = quadraticObjective_->getVectorLengths();
     760  const double *quadraticElement = quadraticObjective_->getElements();
     761  double a = 0.0;
     762  double b = delta;
     763  double c = 0.0;
     764  if (!scaling) {
     765    if (!fullMatrix_) {
     766      int iColumn;
     767      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     768        double valueI = solution[iColumn];
     769        double changeI = change[iColumn];
     770        CoinBigIndex j;
     771        for (j = columnQuadraticStart[iColumn];
     772             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     773          int jColumn = columnQuadratic[j];
     774          double valueJ = solution[jColumn];
     775          double changeJ = change[jColumn];
     776          double elementValue = quadraticElement[j];
     777          if (iColumn != jColumn) {
     778            a += changeI * changeJ * elementValue;
     779            b += (changeI * valueJ + changeJ * valueI) * elementValue;
     780            c += valueI * valueJ * elementValue;
    757781          } else {
    758             COIN_DETAIL_PRINT(printf("odd linear direction %g\n", delta));
    759                return 0.0;
    760           }
    761      }
    762      assert (model);
    763      bool scaling = false;
    764      if ((model->rowScale() ||
    765                model->objectiveScale() != 1.0 || model->optimizationDirection() != 1.0) && inSolve)
    766           scaling = true;
    767      const int * columnQuadratic = quadraticObjective_->getIndices();
    768      const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
    769      const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
    770      const double * quadraticElement = quadraticObjective_->getElements();
    771      double a = 0.0;
    772      double b = delta;
    773      double c = 0.0;
    774      if (!scaling) {
    775           if (!fullMatrix_) {
    776                int iColumn;
    777                for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    778                     double valueI = solution[iColumn];
    779                     double changeI = change[iColumn];
    780                     CoinBigIndex j;
    781                     for (j = columnQuadraticStart[iColumn];
    782                               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    783                          int jColumn = columnQuadratic[j];
    784                          double valueJ = solution[jColumn];
    785                          double changeJ = change[jColumn];
    786                          double elementValue = quadraticElement[j];
    787                          if (iColumn != jColumn) {
    788                               a += changeI * changeJ * elementValue;
    789                               b += (changeI * valueJ + changeJ * valueI) * elementValue;
    790                               c += valueI * valueJ * elementValue;
    791                          } else {
    792                               a += 0.5 * changeI * changeI * elementValue;
    793                               b += changeI * valueI * elementValue;
    794                               c += 0.5 * valueI * valueI * elementValue;
    795                          }
    796                     }
    797                }
     782            a += 0.5 * changeI * changeI * elementValue;
     783            b += changeI * valueI * elementValue;
     784            c += 0.5 * valueI * valueI * elementValue;
     785          }
     786        }
     787      }
     788    } else {
     789      // full matrix stored
     790      int iColumn;
     791      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     792        double valueI = solution[iColumn];
     793        double changeI = change[iColumn];
     794        CoinBigIndex j;
     795        for (j = columnQuadraticStart[iColumn];
     796             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     797          int jColumn = columnQuadratic[j];
     798          double valueJ = solution[jColumn];
     799          double changeJ = change[jColumn];
     800          double elementValue = quadraticElement[j];
     801          valueJ *= elementValue;
     802          a += changeI * changeJ * elementValue;
     803          b += changeI * valueJ;
     804          c += valueI * valueJ;
     805        }
     806      }
     807      a *= 0.5;
     808      c *= 0.5;
     809    }
     810  } else {
     811    // scaling
     812    // for now only if half
     813    assert(!fullMatrix_);
     814    const double *columnScale = model->columnScale();
     815    double direction = model->optimizationDirection() * model->objectiveScale();
     816    // direction is actually scale out not scale in
     817    if (direction)
     818      direction = 1.0 / direction;
     819    if (!columnScale) {
     820      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
     821        double valueI = solution[iColumn];
     822        double changeI = change[iColumn];
     823        CoinBigIndex j;
     824        for (j = columnQuadraticStart[iColumn];
     825             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     826          int jColumn = columnQuadratic[j];
     827          double valueJ = solution[jColumn];
     828          double changeJ = change[jColumn];
     829          double elementValue = quadraticElement[j];
     830          elementValue *= direction;
     831          if (iColumn != jColumn) {
     832            a += changeI * changeJ * elementValue;
     833            b += (changeI * valueJ + changeJ * valueI) * elementValue;
     834            c += valueI * valueJ * elementValue;
    798835          } else {
    799                // full matrix stored
    800                int iColumn;
    801                for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    802                     double valueI = solution[iColumn];
    803                     double changeI = change[iColumn];
    804                     CoinBigIndex j;
    805                     for (j = columnQuadraticStart[iColumn];
    806                               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    807                          int jColumn = columnQuadratic[j];
    808                          double valueJ = solution[jColumn];
    809                          double changeJ = change[jColumn];
    810                          double elementValue = quadraticElement[j];
    811                          valueJ *= elementValue;
    812                          a += changeI * changeJ * elementValue;
    813                          b += changeI * valueJ;
    814                          c += valueI * valueJ;
    815                     }
    816                }
    817                a *= 0.5;
    818                c *= 0.5;
    819           }
    820      } else {
    821           // scaling
    822           // for now only if half
    823           assert (!fullMatrix_);
    824           const double * columnScale = model->columnScale();
    825           double direction = model->optimizationDirection() * model->objectiveScale();
    826           // direction is actually scale out not scale in
    827           if (direction)
    828                direction = 1.0 / direction;
    829           if (!columnScale) {
    830                for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
    831                     double valueI = solution[iColumn];
    832                     double changeI = change[iColumn];
    833                     CoinBigIndex j;
    834                     for (j = columnQuadraticStart[iColumn];
    835                               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    836                          int jColumn = columnQuadratic[j];
    837                          double valueJ = solution[jColumn];
    838                          double changeJ = change[jColumn];
    839                          double elementValue = quadraticElement[j];
    840                          elementValue *= direction;
    841                          if (iColumn != jColumn) {
    842                               a += changeI * changeJ * elementValue;
    843                               b += (changeI * valueJ + changeJ * valueI) * elementValue;
    844                               c += valueI * valueJ * elementValue;
    845                          } else {
    846                               a += 0.5 * changeI * changeI * elementValue;
    847                               b += changeI * valueI * elementValue;
    848                               c += 0.5 * valueI * valueI * elementValue;
    849                          }
    850                     }
    851                }
     836            a += 0.5 * changeI * changeI * elementValue;
     837            b += changeI * valueI * elementValue;
     838            c += 0.5 * valueI * valueI * elementValue;
     839          }
     840        }
     841      }
     842    } else {
     843      // scaling
     844      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
     845        double valueI = solution[iColumn];
     846        double changeI = change[iColumn];
     847        double scaleI = columnScale[iColumn] * direction;
     848        CoinBigIndex j;
     849        for (j = columnQuadraticStart[iColumn];
     850             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     851          int jColumn = columnQuadratic[j];
     852          double valueJ = solution[jColumn];
     853          double changeJ = change[jColumn];
     854          double elementValue = quadraticElement[j];
     855          elementValue *= scaleI * columnScale[jColumn];
     856          if (iColumn != jColumn) {
     857            a += changeI * changeJ * elementValue;
     858            b += (changeI * valueJ + changeJ * valueI) * elementValue;
     859            c += valueI * valueJ * elementValue;
    852860          } else {
    853                // scaling
    854                for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
    855                     double valueI = solution[iColumn];
    856                     double changeI = change[iColumn];
    857                     double scaleI = columnScale[iColumn] * direction;
    858                     CoinBigIndex j;
    859                     for (j = columnQuadraticStart[iColumn];
    860                               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    861                          int jColumn = columnQuadratic[j];
    862                          double valueJ = solution[jColumn];
    863                          double changeJ = change[jColumn];
    864                          double elementValue = quadraticElement[j];
    865                          elementValue *= scaleI * columnScale[jColumn];
    866                          if (iColumn != jColumn) {
    867                               a += changeI * changeJ * elementValue;
    868                               b += (changeI * valueJ + changeJ * valueI) * elementValue;
    869                               c += valueI * valueJ * elementValue;
    870                          } else {
    871                               a += 0.5 * changeI * changeI * elementValue;
    872                               b += changeI * valueI * elementValue;
    873                               c += 0.5 * valueI * valueI * elementValue;
    874                          }
    875                     }
    876                }
    877           }
    878      }
    879      double theta;
    880      //printf("Current cost %g\n",c+linearCost);
    881      currentObj = c + linearCost;
    882      thetaObj = currentObj + a * maximumTheta * maximumTheta + b * maximumTheta;
    883      // minimize a*x*x + b*x + c
    884      if (a <= 0.0) {
    885           theta = maximumTheta;
    886      } else {
    887           theta = -0.5 * b / a;
    888      }
    889      predictedObj = currentObj + a * theta * theta + b * theta;
    890      if (b > 0.0) {
    891           if (model->messageHandler()->logLevel() & 32)
    892                printf("a %g b %g c %g => %g\n", a, b, c, theta);
    893           b = 0.0;
    894      }
    895      return CoinMin(theta, maximumTheta);
     861            a += 0.5 * changeI * changeI * elementValue;
     862            b += changeI * valueI * elementValue;
     863            c += 0.5 * valueI * valueI * elementValue;
     864          }
     865        }
     866      }
     867    }
     868  }
     869  double theta;
     870  //printf("Current cost %g\n",c+linearCost);
     871  currentObj = c + linearCost;
     872  thetaObj = currentObj + a * maximumTheta * maximumTheta + b * maximumTheta;
     873  // minimize a*x*x + b*x + c
     874  if (a <= 0.0) {
     875    theta = maximumTheta;
     876  } else {
     877    theta = -0.5 * b / a;
     878  }
     879  predictedObj = currentObj + a * theta * theta + b * theta;
     880  if (b > 0.0) {
     881    if (model->messageHandler()->logLevel() & 32)
     882      printf("a %g b %g c %g => %g\n", a, b, c, theta);
     883    b = 0.0;
     884  }
     885  return CoinMin(theta, maximumTheta);
    896886}
    897887// Return objective value (without any ClpModel offset) (model may be NULL)
    898888double
    899 ClpQuadraticObjective::objectiveValue(const ClpSimplex * model, const double * solution) const
    900 {
    901      bool scaling = false;
    902      if (model && (model->rowScale() ||
    903                    model->objectiveScale() != 1.0))
    904           scaling = true;
    905      const double * cost = NULL;
    906      if (model)
    907           cost = model->costRegion();
    908      if (!cost) {
    909           // not in solve
    910           cost = objective_;
    911           scaling = false;
    912      }
    913      double linearCost = 0.0;
    914      int numberColumns = model->numberColumns();
    915      int numberTotal = numberColumns;
    916      double currentObj = 0.0;
    917      for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
    918           linearCost += cost[iColumn] * solution[iColumn];
    919      }
    920      if (!activated_ || !quadraticObjective_) {
    921           currentObj = linearCost;
    922           return currentObj;
    923      }
    924      assert (model);
    925      const int * columnQuadratic = quadraticObjective_->getIndices();
    926      const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
    927      const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
    928      const double * quadraticElement = quadraticObjective_->getElements();
    929      double c = 0.0;
    930      if (!scaling) {
    931           if (!fullMatrix_) {
    932                int iColumn;
    933                for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    934                     double valueI = solution[iColumn];
    935                     CoinBigIndex j;
    936                     for (j = columnQuadraticStart[iColumn];
    937                               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    938                          int jColumn = columnQuadratic[j];
    939                          double valueJ = solution[jColumn];
    940                          double elementValue = quadraticElement[j];
    941                          if (iColumn != jColumn) {
    942                               c += valueI * valueJ * elementValue;
    943                          } else {
    944                               c += 0.5 * valueI * valueI * elementValue;
    945                          }
    946                     }
    947                }
     889ClpQuadraticObjective::objectiveValue(const ClpSimplex *model, const double *solution) const
     890{
     891  bool scaling = false;
     892  if (model && (model->rowScale() || model->objectiveScale() != 1.0))
     893    scaling = true;
     894  const double *cost = NULL;
     895  if (model)
     896    cost = model->costRegion();
     897  if (!cost) {
     898    // not in solve
     899    cost = objective_;
     900    scaling = false;
     901  }
     902  double linearCost = 0.0;
     903  int numberColumns = model->numberColumns();
     904  int numberTotal = numberColumns;
     905  double currentObj = 0.0;
     906  for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
     907    linearCost += cost[iColumn] * solution[iColumn];
     908  }
     909  if (!activated_ || !quadraticObjective_) {
     910    currentObj = linearCost;
     911    return currentObj;
     912  }
     913  assert(model);
     914  const int *columnQuadratic = quadraticObjective_->getIndices();
     915  const CoinBigIndex *columnQuadraticStart = quadraticObjective_->getVectorStarts();
     916  const int *columnQuadraticLength = quadraticObjective_->getVectorLengths();
     917  const double *quadraticElement = quadraticObjective_->getElements();
     918  double c = 0.0;
     919  if (!scaling) {
     920    if (!fullMatrix_) {
     921      int iColumn;
     922      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     923        double valueI = solution[iColumn];
     924        CoinBigIndex j;
     925        for (j = columnQuadraticStart[iColumn];
     926             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     927          int jColumn = columnQuadratic[j];
     928          double valueJ = solution[jColumn];
     929          double elementValue = quadraticElement[j];
     930          if (iColumn != jColumn) {
     931            c += valueI * valueJ * elementValue;
    948932          } else {
    949                // full matrix stored
    950                int iColumn;
    951                for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    952                     double valueI = solution[iColumn];
    953                     CoinBigIndex j;
    954                     for (j = columnQuadraticStart[iColumn];
    955                               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    956                          int jColumn = columnQuadratic[j];
    957                          double valueJ = solution[jColumn];
    958                          double elementValue = quadraticElement[j];
    959                          valueJ *= elementValue;
    960                          c += valueI * valueJ;
    961                     }
    962                }
    963                c *= 0.5;
    964           }
    965      } else {
    966           // scaling
    967           // for now only if half
    968           assert (!fullMatrix_);
    969           const double * columnScale = model->columnScale();
    970           double direction = model->objectiveScale();
    971           // direction is actually scale out not scale in
    972           if (direction)
    973                direction = 1.0 / direction;
    974           if (!columnScale) {
    975                for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
    976                     double valueI = solution[iColumn];
    977                     CoinBigIndex j;
    978                     for (j = columnQuadraticStart[iColumn];
    979                               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    980                          int jColumn = columnQuadratic[j];
    981                          double valueJ = solution[jColumn];
    982                          double elementValue = quadraticElement[j];
    983                          elementValue *= direction;
    984                          if (iColumn != jColumn) {
    985                               c += valueI * valueJ * elementValue;
    986                          } else {
    987                               c += 0.5 * valueI * valueI * elementValue;
    988                          }
    989                     }
    990                }
     933            c += 0.5 * valueI * valueI * elementValue;
     934          }
     935        }
     936      }
     937    } else {
     938      // full matrix stored
     939      int iColumn;
     940      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     941        double valueI = solution[iColumn];
     942        CoinBigIndex j;
     943        for (j = columnQuadraticStart[iColumn];
     944             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     945          int jColumn = columnQuadratic[j];
     946          double valueJ = solution[jColumn];
     947          double elementValue = quadraticElement[j];
     948          valueJ *= elementValue;
     949          c += valueI * valueJ;
     950        }
     951      }
     952      c *= 0.5;
     953    }
     954  } else {
     955    // scaling
     956    // for now only if half
     957    assert(!fullMatrix_);
     958    const double *columnScale = model->columnScale();
     959    double direction = model->objectiveScale();
     960    // direction is actually scale out not scale in
     961    if (direction)
     962      direction = 1.0 / direction;
     963    if (!columnScale) {
     964      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
     965        double valueI = solution[iColumn];
     966        CoinBigIndex j;
     967        for (j = columnQuadraticStart[iColumn];
     968             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     969          int jColumn = columnQuadratic[j];
     970          double valueJ = solution[jColumn];
     971          double elementValue = quadraticElement[j];
     972          elementValue *= direction;
     973          if (iColumn != jColumn) {
     974            c += valueI * valueJ * elementValue;
    991975          } else {
    992                // scaling
    993                for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
    994                     double valueI = solution[iColumn];
    995                     double scaleI = columnScale[iColumn] * direction;
    996                     CoinBigIndex j;
    997                     for (j = columnQuadraticStart[iColumn];
    998                               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    999                          int jColumn = columnQuadratic[j];
    1000                          double valueJ = solution[jColumn];
    1001                          double elementValue = quadraticElement[j];
    1002                          elementValue *= scaleI * columnScale[jColumn];
    1003                          if (iColumn != jColumn) {
    1004                               c += valueI * valueJ * elementValue;
    1005                          } else {
    1006                               c += 0.5 * valueI * valueI * elementValue;
    1007                          }
    1008                     }
    1009                }
    1010           }
    1011      }
    1012      currentObj = c + linearCost;
    1013      return currentObj;
     976            c += 0.5 * valueI * valueI * elementValue;
     977          }
     978        }
     979      }
     980    } else {
     981      // scaling
     982      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
     983        double valueI = solution[iColumn];
     984        double scaleI = columnScale[iColumn] * direction;
     985        CoinBigIndex j;
     986        for (j = columnQuadraticStart[iColumn];
     987             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     988          int jColumn = columnQuadratic[j];
     989          double valueJ = solution[jColumn];
     990          double elementValue = quadraticElement[j];
     991          elementValue *= scaleI * columnScale[jColumn];
     992          if (iColumn != jColumn) {
     993            c += valueI * valueJ * elementValue;
     994          } else {
     995            c += 0.5 * valueI * valueI * elementValue;
     996          }
     997        }
     998      }
     999    }
     1000  }
     1001  currentObj = c + linearCost;
     1002  return currentObj;
    10141003}
    10151004// Scale objective
    1016 void
    1017 ClpQuadraticObjective::reallyScale(const double * columnScale)
    1018 {
    1019      const int * columnQuadratic = quadraticObjective_->getIndices();
    1020      const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
    1021      const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
    1022      double * quadraticElement = quadraticObjective_->getMutableElements();
    1023      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1024           double scaleI = columnScale[iColumn];
    1025           objective_[iColumn] *= scaleI;
    1026           CoinBigIndex j;
    1027           for (j = columnQuadraticStart[iColumn];
    1028                     j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    1029                int jColumn = columnQuadratic[j];
    1030                quadraticElement[j] *= scaleI * columnScale[jColumn];
    1031           }
    1032      }
     1005void ClpQuadraticObjective::reallyScale(const double *columnScale)
     1006{
     1007  const int *columnQuadratic = quadraticObjective_->getIndices();
     1008  const CoinBigIndex *columnQuadraticStart = quadraticObjective_->getVectorStarts();
     1009  const int *columnQuadraticLength = quadraticObjective_->getVectorLengths();
     1010  double *quadraticElement = quadraticObjective_->getMutableElements();
     1011  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
     1012    double scaleI = columnScale[iColumn];
     1013    objective_[iColumn] *= scaleI;
     1014    CoinBigIndex j;
     1015    for (j = columnQuadraticStart[iColumn];
     1016         j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     1017      int jColumn = columnQuadratic[j];
     1018      quadraticElement[j] *= scaleI * columnScale[jColumn];
     1019    }
     1020  }
    10331021}
    10341022/* Given a zeroed array sets nonlinear columns to 1.
    10351023   Returns number of nonlinear columns
    10361024*/
    1037 int
    1038 ClpQuadraticObjective::markNonlinear(char * which)
    1039 {
    1040      int iColumn;
    1041      const int * columnQuadratic = quadraticObjective_->getIndices();
    1042      const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
    1043      const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
    1044      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1045           CoinBigIndex j;
    1046           for (j = columnQuadraticStart[iColumn];
    1047                     j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    1048                int jColumn = columnQuadratic[j];
    1049                which[jColumn] = 1;
    1050                which[iColumn] = 1;
    1051           }
    1052      }
    1053      int numberNonLinearColumns = 0;
    1054      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1055           if(which[iColumn])
    1056                numberNonLinearColumns++;
    1057      }
    1058      return numberNonLinearColumns;
    1059 }
     1025int ClpQuadraticObjective::markNonlinear(char *which)
     1026{
     1027  int iColumn;
     1028  const int *columnQuadratic = quadraticObjective_->getIndices();
     1029  const CoinBigIndex *columnQuadraticStart = quadraticObjective_->getVectorStarts();
     1030  const int *columnQuadraticLength = quadraticObjective_->getVectorLengths();
     1031  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     1032    CoinBigIndex j;
     1033    for (j = columnQuadraticStart[iColumn];
     1034         j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     1035      int jColumn = columnQuadratic[j];
     1036      which[jColumn] = 1;
     1037      which[iColumn] = 1;
     1038    }
     1039  }
     1040  int numberNonLinearColumns = 0;
     1041  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     1042    if (which[iColumn])
     1043      numberNonLinearColumns++;
     1044  }
     1045  return numberNonLinearColumns;
     1046}
     1047
     1048/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     1049*/
Note: See TracChangeset for help on using the changeset viewer.