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

formatting

File:
1 edited

Legend:

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

    r1723 r2385  
    44// This code is licensed under the terms of the Eclipse Public License (EPL).
    55
    6 
    76#include "CoinPragma.hpp"
    87#include "CoinHelperFunctions.hpp"
     
    2120// Default Constructor
    2221//-------------------------------------------------------------------
    23 ClpCholeskyWssmpKKT::ClpCholeskyWssmpKKT (int denseThreshold)
    24      : ClpCholeskyBase(denseThreshold)
    25 {
    26      type_ = 21;
     22ClpCholeskyWssmpKKT::ClpCholeskyWssmpKKT(int denseThreshold)
     23  : ClpCholeskyBase(denseThreshold)
     24{
     25  type_ = 21;
    2726}
    2827
     
    3029// Copy constructor
    3130//-------------------------------------------------------------------
    32 ClpCholeskyWssmpKKT::ClpCholeskyWssmpKKT (const ClpCholeskyWssmpKKT & rhs)
    33      : ClpCholeskyBase(rhs)
    34 {
    35 }
    36 
     31ClpCholeskyWssmpKKT::ClpCholeskyWssmpKKT(const ClpCholeskyWssmpKKT &rhs)
     32  : ClpCholeskyBase(rhs)
     33{
     34}
    3735
    3836//-------------------------------------------------------------------
    3937// Destructor
    4038//-------------------------------------------------------------------
    41 ClpCholeskyWssmpKKT::~ClpCholeskyWssmpKKT ()
     39ClpCholeskyWssmpKKT::~ClpCholeskyWssmpKKT()
    4240{
    4341}
     
    4745//-------------------------------------------------------------------
    4846ClpCholeskyWssmpKKT &
    49 ClpCholeskyWssmpKKT::operator=(const ClpCholeskyWssmpKKT& rhs)
    50 {
    51      if (this != &rhs) {
    52           ClpCholeskyBase::operator=(rhs);
    53      }
    54      return *this;
     47ClpCholeskyWssmpKKT::operator=(const ClpCholeskyWssmpKKT &rhs)
     48{
     49  if (this != &rhs) {
     50    ClpCholeskyBase::operator=(rhs);
     51  }
     52  return *this;
    5553}
    5654//-------------------------------------------------------------------
    5755// Clone
    5856//-------------------------------------------------------------------
    59 ClpCholeskyBase * ClpCholeskyWssmpKKT::clone() const
    60 {
    61      return new ClpCholeskyWssmpKKT(*this);
     57ClpCholeskyBase *ClpCholeskyWssmpKKT::clone() const
     58{
     59  return new ClpCholeskyWssmpKKT(*this);
    6260}
    6361// At present I can't get wssmp to work as my libraries seem to be out of sync
     
    6563#ifndef USE_EKKWSSMP
    6664extern "C" {
    67 void F77_FUNC(wsetmaxthrds,WSETMAXTHRDS)(const int* NTHREADS);
    68 
    69 void F77_FUNC(wssmp,WSSMP)(const int* N, const int* IA,
    70                            const int* JA, const double* AVALS,
    71                            double* DIAG,  int* PERM,
    72                            int* INVP,  double* B,   
    73                            const int* LDB, const int* NRHS,
    74                            double* AUX, const int* NAUX,
    75                            int* MRP, int* IPARM,
    76                            double* DPARM);
    77 void F77_FUNC_(wsmp_clear,WSMP_CLEAR)(void);
     65void F77_FUNC(wsetmaxthrds, WSETMAXTHRDS)(const int *NTHREADS);
     66
     67void F77_FUNC(wssmp, WSSMP)(const int *N, const int *IA,
     68  const int *JA, const double *AVALS,
     69  double *DIAG, int *PERM,
     70  int *INVP, double *B,
     71  const int *LDB, const int *NRHS,
     72  double *AUX, const int *NAUX,
     73  int *MRP, int *IPARM,
     74  double *DPARM);
     75void F77_FUNC_(wsmp_clear, WSMP_CLEAR)(void);
    7876}
    7977#else
     
    8280typedef struct EKKContext EKKContext;
    8381
    84 
    8582extern "C" {
    86      EKKContext *  ekk_initializeContext();
    87      void ekk_endContext(EKKContext * context);
    88      EKKModel *  ekk_newModel(EKKContext * env, const char * name);
    89      int ekk_deleteModel(EKKModel * model);
    90 }
    91 static  EKKModel * model = NULL;
    92 static  EKKContext * context = NULL;
    93 extern "C" void ekkwssmp(EKKModel *, int * n,
    94                          int * columnStart , int * rowIndex , double * element,
    95                          double * diagonal , int * perm , int * invp ,
    96                          double * rhs , int * ldb , int * nrhs ,
    97                          double * aux , int * naux ,
    98                          int   * mrp , int * iparm , double * dparm);
    99 static void F77_FUNC(wssmp,WSSMP)( int *n, int *ia, int *ja,
    100                    double *avals, double *diag, int *perm, int *invp,
    101                    double *b, int *ldb, int *nrhs, double *aux, int *
    102                    naux, int *mrp, int *iparm, double *dparm)
    103 {
    104      if (!context) {
    105           /* initialize OSL environment */
    106           context = ekk_initializeContext();
    107           model = ekk_newModel(context, "");
    108      }
    109      ekkwssmp(model, n, ia, ja,
    110               avals, diag, perm, invp,
    111               b, ldb, nrhs, aux,
    112               naux, mrp, iparm, dparm);
    113      //ekk_deleteModel(model);
    114      //ekk_endContext(context);
     83EKKContext *ekk_initializeContext();
     84void ekk_endContext(EKKContext *context);
     85EKKModel *ekk_newModel(EKKContext *env, const char *name);
     86int ekk_deleteModel(EKKModel *model);
     87}
     88static EKKModel *model = NULL;
     89static EKKContext *context = NULL;
     90extern "C" void ekkwssmp(EKKModel *, int *n,
     91  int *columnStart, int *rowIndex, double *element,
     92  double *diagonal, int *perm, int *invp,
     93  double *rhs, int *ldb, int *nrhs,
     94  double *aux, int *naux,
     95  int *mrp, int *iparm, double *dparm);
     96static void F77_FUNC(wssmp, WSSMP)(int *n, int *ia, int *ja,
     97  double *avals, double *diag, int *perm, int *invp,
     98  double *b, int *ldb, int *nrhs, double *aux, int *naux, int *mrp, int *iparm, double *dparm)
     99{
     100  if (!context) {
     101    /* initialize OSL environment */
     102    context = ekk_initializeContext();
     103    model = ekk_newModel(context, "");
     104  }
     105  ekkwssmp(model, n, ia, ja,
     106    avals, diag, perm, invp,
     107    b, ldb, nrhs, aux,
     108    naux, mrp, iparm, dparm);
     109  //ekk_deleteModel(model);
     110  //ekk_endContext(context);
    115111}
    116112#endif
    117113
    118114/* Orders rows and saves pointer to model */
    119 int
    120 ClpCholeskyWssmpKKT::order(ClpInterior * model)
    121 {
    122      int numberRowsModel = model->numberRows();
    123      int numberColumns = model->numberColumns();
    124      int numberTotal = numberColumns + numberRowsModel;
    125      numberRows_ = 2 * numberRowsModel + numberColumns;
    126      rowsDropped_ = new char [numberRows_];
    127      memset(rowsDropped_, 0, numberRows_);
    128      numberRowsDropped_ = 0;
    129      model_ = model;
    130      CoinPackedMatrix * quadratic = NULL;
    131      ClpQuadraticObjective * quadraticObj =
    132           (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
    133      if (quadraticObj)
    134           quadratic = quadraticObj->quadraticObjective();
    135      int numberElements = model_->clpMatrix()->getNumElements();
    136      numberElements = numberElements + 2 * numberRowsModel + numberTotal;
    137      if (quadratic)
    138           numberElements += quadratic->getNumElements();
    139      // Space for starts
    140      choleskyStart_ = new CoinBigIndex[numberRows_+1];
    141      const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    142      const int * columnLength = model_->clpMatrix()->getVectorLengths();
    143      const int * row = model_->clpMatrix()->getIndices();
    144      //const double * element = model_->clpMatrix()->getElements();
    145      // Now we have size - create arrays and fill in
    146      try {
    147           choleskyRow_ = new int [numberElements];
    148      } catch (...) {
    149           // no memory
    150           delete [] choleskyStart_;
    151           choleskyStart_ = NULL;
    152           return -1;
    153      }
    154      try {
    155           sparseFactor_ = new double[numberElements];
    156      } catch (...) {
    157           // no memory
    158           delete [] choleskyRow_;
    159           choleskyRow_ = NULL;
    160           delete [] choleskyStart_;
    161           choleskyStart_ = NULL;
    162           return -1;
    163      }
    164      int iRow, iColumn;
    165 
    166      sizeFactor_ = 0;
    167      // matrix
    168      if (!quadratic) {
    169           for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    170                choleskyStart_[iColumn] = sizeFactor_;
    171                choleskyRow_[sizeFactor_++] = iColumn;
    172                CoinBigIndex start = columnStart[iColumn];
    173                CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    174                for (CoinBigIndex j = start; j < end; j++) {
    175                     choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
    176                }
    177           }
    178      } else {
    179           // Quadratic
    180           const int * columnQuadratic = quadratic->getIndices();
    181           const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
    182           const int * columnQuadraticLength = quadratic->getVectorLengths();
    183           //const double * quadraticElement = quadratic->getElements();
    184           for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    185                choleskyStart_[iColumn] = sizeFactor_;
    186                choleskyRow_[sizeFactor_++] = iColumn;
    187                for (CoinBigIndex j = columnQuadraticStart[iColumn];
    188                          j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    189                     int jColumn = columnQuadratic[j];
    190                     if (jColumn > iColumn)
    191                          choleskyRow_[sizeFactor_++] = jColumn;
    192                }
    193                CoinBigIndex start = columnStart[iColumn];
    194                CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    195                for (CoinBigIndex j = start; j < end; j++) {
    196                     choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
    197                }
    198           }
    199      }
    200      // slacks
    201      for (; iColumn < numberTotal; iColumn++) {
    202           choleskyStart_[iColumn] = sizeFactor_;
    203           choleskyRow_[sizeFactor_++] = iColumn;
    204           choleskyRow_[sizeFactor_++] = iColumn - numberColumns + numberTotal;
    205      }
    206      // Transpose - nonzero diagonal (may regularize)
    207      for (iRow = 0; iRow < numberRowsModel; iRow++) {
    208           choleskyStart_[iRow+numberTotal] = sizeFactor_;
    209           // diagonal
    210           choleskyRow_[sizeFactor_++] = iRow + numberTotal;
    211      }
    212      choleskyStart_[numberRows_] = sizeFactor_;
    213      permuteInverse_ = new int [numberRows_];
    214      permute_ = new int[numberRows_];
    215      integerParameters_[0] = 0;
    216      int i0 = 0;
    217      int i1 = 1;
     115int ClpCholeskyWssmpKKT::order(ClpInterior *model)
     116{
     117  int numberRowsModel = model->numberRows();
     118  int numberColumns = model->numberColumns();
     119  int numberTotal = numberColumns + numberRowsModel;
     120  numberRows_ = 2 * numberRowsModel + numberColumns;
     121  rowsDropped_ = new char[numberRows_];
     122  memset(rowsDropped_, 0, numberRows_);
     123  numberRowsDropped_ = 0;
     124  model_ = model;
     125  CoinPackedMatrix *quadratic = NULL;
     126  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(model_->objectiveAsObject()));
     127  if (quadraticObj)
     128    quadratic = quadraticObj->quadraticObjective();
     129  int numberElements = model_->clpMatrix()->getNumElements();
     130  numberElements = numberElements + 2 * numberRowsModel + numberTotal;
     131  if (quadratic)
     132    numberElements += quadratic->getNumElements();
     133  // Space for starts
     134  choleskyStart_ = new CoinBigIndex[numberRows_ + 1];
     135  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
     136  const int *columnLength = model_->clpMatrix()->getVectorLengths();
     137  const int *row = model_->clpMatrix()->getIndices();
     138  //const double * element = model_->clpMatrix()->getElements();
     139  // Now we have size - create arrays and fill in
     140  try {
     141    choleskyRow_ = new int[numberElements];
     142  } catch (...) {
     143    // no memory
     144    delete[] choleskyStart_;
     145    choleskyStart_ = NULL;
     146    return -1;
     147  }
     148  try {
     149    sparseFactor_ = new double[numberElements];
     150  } catch (...) {
     151    // no memory
     152    delete[] choleskyRow_;
     153    choleskyRow_ = NULL;
     154    delete[] choleskyStart_;
     155    choleskyStart_ = NULL;
     156    return -1;
     157  }
     158  int iRow, iColumn;
     159
     160  sizeFactor_ = 0;
     161  // matrix
     162  if (!quadratic) {
     163    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     164      choleskyStart_[iColumn] = sizeFactor_;
     165      choleskyRow_[sizeFactor_++] = iColumn;
     166      CoinBigIndex start = columnStart[iColumn];
     167      CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     168      for (CoinBigIndex j = start; j < end; j++) {
     169        choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
     170      }
     171    }
     172  } else {
     173    // Quadratic
     174    const int *columnQuadratic = quadratic->getIndices();
     175    const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
     176    const int *columnQuadraticLength = quadratic->getVectorLengths();
     177    //const double * quadraticElement = quadratic->getElements();
     178    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     179      choleskyStart_[iColumn] = sizeFactor_;
     180      choleskyRow_[sizeFactor_++] = iColumn;
     181      for (CoinBigIndex j = columnQuadraticStart[iColumn];
     182           j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     183        int jColumn = columnQuadratic[j];
     184        if (jColumn > iColumn)
     185          choleskyRow_[sizeFactor_++] = jColumn;
     186      }
     187      CoinBigIndex start = columnStart[iColumn];
     188      CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     189      for (CoinBigIndex j = start; j < end; j++) {
     190        choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
     191      }
     192    }
     193  }
     194  // slacks
     195  for (; iColumn < numberTotal; iColumn++) {
     196    choleskyStart_[iColumn] = sizeFactor_;
     197    choleskyRow_[sizeFactor_++] = iColumn;
     198    choleskyRow_[sizeFactor_++] = iColumn - numberColumns + numberTotal;
     199  }
     200  // Transpose - nonzero diagonal (may regularize)
     201  for (iRow = 0; iRow < numberRowsModel; iRow++) {
     202    choleskyStart_[iRow + numberTotal] = sizeFactor_;
     203    // diagonal
     204    choleskyRow_[sizeFactor_++] = iRow + numberTotal;
     205  }
     206  choleskyStart_[numberRows_] = sizeFactor_;
     207  permuteInverse_ = new int[numberRows_];
     208  permute_ = new int[numberRows_];
     209  integerParameters_[0] = 0;
     210  int i0 = 0;
     211  int i1 = 1;
    218212#ifndef USE_EKKWSSMP
    219      int i2 = 1;
    220      if (model->numberThreads() <= 0)
    221           i2 = 1;
    222      else
    223           i2 = model->numberThreads();
    224      F77_FUNC(wsetmaxthrds,WSETMAXTHRDS)(&i2);
    225 #endif
    226      F77_FUNC(wssmp,WSSMP)(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
    227            NULL, permute_, permuteInverse_, 0, &numberRows_, &i1,
    228            NULL, &i0, NULL, integerParameters_, doubleParameters_);
    229      integerParameters_[1] = 1; //order and symbolic
    230      integerParameters_[2] = 2;
    231      integerParameters_[3] = 0; //CSR
    232      integerParameters_[4] = 0; //C style
    233      integerParameters_[13] = 1; //reuse initial factorization space
    234      integerParameters_[15+0] = 1; //ordering
    235      integerParameters_[15+1] = 0;
    236      integerParameters_[15+2] = 1;
    237      integerParameters_[15+3] = 0;
    238      integerParameters_[15+4] = 1;
    239      doubleParameters_[10] = 1.0e-20;
    240      doubleParameters_[11] = 1.0e-15;
     213  int i2 = 1;
     214  if (model->numberThreads() <= 0)
     215    i2 = 1;
     216  else
     217    i2 = model->numberThreads();
     218  F77_FUNC(wsetmaxthrds, WSETMAXTHRDS)
     219  (&i2);
     220#endif
     221  F77_FUNC(wssmp, WSSMP)
     222  (&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
     223    NULL, permute_, permuteInverse_, 0, &numberRows_, &i1,
     224    NULL, &i0, NULL, integerParameters_, doubleParameters_);
     225  integerParameters_[1] = 1; //order and symbolic
     226  integerParameters_[2] = 2;
     227  integerParameters_[3] = 0; //CSR
     228  integerParameters_[4] = 0; //C style
     229  integerParameters_[13] = 1; //reuse initial factorization space
     230  integerParameters_[15 + 0] = 1; //ordering
     231  integerParameters_[15 + 1] = 0;
     232  integerParameters_[15 + 2] = 1;
     233  integerParameters_[15 + 3] = 0;
     234  integerParameters_[15 + 4] = 1;
     235  doubleParameters_[10] = 1.0e-20;
     236  doubleParameters_[11] = 1.0e-15;
    241237#if 1
    242      integerParameters_[1] = 2; //just symbolic
    243      for (int iRow = 0; iRow < numberRows_; iRow++) {
    244           permuteInverse_[iRow] = iRow;
    245           permute_[iRow] = iRow;
    246      }
    247 #endif
    248      F77_FUNC(wssmp,WSSMP)(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
    249            NULL, permute_, permuteInverse_, NULL, &numberRows_, &i1,
    250            NULL, &i0, NULL, integerParameters_, doubleParameters_);
    251      //std::cout<<"Ordering and symbolic factorization took "<<doubleParameters_[0]<<std::endl;
    252      if (integerParameters_[63]) {
    253           std::cout << "wssmp returning error code of " << integerParameters_[63] << std::endl;
    254           return 1;
    255      }
    256      std::cout << integerParameters_[23] << " elements in sparse Cholesky" << std::endl;
    257      if (!integerParameters_[23]) {
    258           for (int iRow = 0; iRow < numberRows_; iRow++) {
    259                permuteInverse_[iRow] = iRow;
    260                permute_[iRow] = iRow;
    261           }
    262           std::cout << "wssmp says no elements - fully dense? - switching to dense" << std::endl;
    263           integerParameters_[1] = 2;
    264           integerParameters_[2] = 2;
    265           integerParameters_[7] = 1; // no permute
    266           F77_FUNC(wssmp,WSSMP)(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
    267                 NULL, permute_, permuteInverse_, NULL, &numberRows_, &i1,
    268                 NULL, &i0, NULL, integerParameters_, doubleParameters_);
    269           std::cout << integerParameters_[23] << " elements in dense Cholesky" << std::endl;
    270      }
    271      return 0;
     238  integerParameters_[1] = 2; //just symbolic
     239  for (int iRow = 0; iRow < numberRows_; iRow++) {
     240    permuteInverse_[iRow] = iRow;
     241    permute_[iRow] = iRow;
     242  }
     243#endif
     244  F77_FUNC(wssmp, WSSMP)
     245  (&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
     246    NULL, permute_, permuteInverse_, NULL, &numberRows_, &i1,
     247    NULL, &i0, NULL, integerParameters_, doubleParameters_);
     248  //std::cout<<"Ordering and symbolic factorization took "<<doubleParameters_[0]<<std::endl;
     249  if (integerParameters_[63]) {
     250    std::cout << "wssmp returning error code of " << integerParameters_[63] << std::endl;
     251    return 1;
     252  }
     253  std::cout << integerParameters_[23] << " elements in sparse Cholesky" << std::endl;
     254  if (!integerParameters_[23]) {
     255    for (int iRow = 0; iRow < numberRows_; iRow++) {
     256      permuteInverse_[iRow] = iRow;
     257      permute_[iRow] = iRow;
     258    }
     259    std::cout << "wssmp says no elements - fully dense? - switching to dense" << std::endl;
     260    integerParameters_[1] = 2;
     261    integerParameters_[2] = 2;
     262    integerParameters_[7] = 1; // no permute
     263    F77_FUNC(wssmp, WSSMP)
     264    (&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
     265      NULL, permute_, permuteInverse_, NULL, &numberRows_, &i1,
     266      NULL, &i0, NULL, integerParameters_, doubleParameters_);
     267    std::cout << integerParameters_[23] << " elements in dense Cholesky" << std::endl;
     268  }
     269  return 0;
    272270}
    273271/* Does Symbolic factorization given permutation.
     
    275273   user must provide factorize and solve.  Otherwise the default factorization is used
    276274   returns non-zero if not enough memory */
    277 int
    278 ClpCholeskyWssmpKKT::symbolic()
    279 {
    280      return 0;
     275int ClpCholeskyWssmpKKT::symbolic()
     276{
     277  return 0;
    281278}
    282279/* Factorize - filling in rowsDropped and returning number dropped */
    283 int
    284 ClpCholeskyWssmpKKT::factorize(const double * diagonal, int * rowsDropped)
    285 {
    286      int numberRowsModel = model_->numberRows();
    287      int numberColumns = model_->numberColumns();
    288      int numberTotal = numberColumns + numberRowsModel;
    289      int newDropped = 0;
    290      double largest = 0.0;
    291      double smallest;
    292      //perturbation
    293      double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
    294      perturbation = perturbation * perturbation;
    295      if (perturbation > 1.0) {
     280int ClpCholeskyWssmpKKT::factorize(const double *diagonal, int *rowsDropped)
     281{
     282  int numberRowsModel = model_->numberRows();
     283  int numberColumns = model_->numberColumns();
     284  int numberTotal = numberColumns + numberRowsModel;
     285  int newDropped = 0;
     286  double largest = 0.0;
     287  double smallest;
     288  //perturbation
     289  double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
     290  perturbation = perturbation * perturbation;
     291  if (perturbation > 1.0) {
    296292#ifdef COIN_DEVELOP
    297           //if (model_->model()->logLevel()&4)
    298           std::cout << "large perturbation " << perturbation << std::endl;
    299 #endif
    300           perturbation = sqrt(perturbation);;
    301           perturbation = 1.0;
    302      }
    303      // need to recreate every time
    304      int iRow, iColumn;
    305      const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    306      const int * columnLength = model_->clpMatrix()->getVectorLengths();
    307      const int * row = model_->clpMatrix()->getIndices();
    308      const double * element = model_->clpMatrix()->getElements();
    309 
    310      CoinBigIndex numberElements = 0;
    311      CoinPackedMatrix * quadratic = NULL;
    312      ClpQuadraticObjective * quadraticObj =
    313           (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
    314      if (quadraticObj)
    315           quadratic = quadraticObj->quadraticObjective();
    316      // matrix
    317      if (!quadratic) {
    318           for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    319                choleskyStart_[iColumn] = numberElements;
    320                double value = diagonal[iColumn];
    321                if (fabs(value) > 1.0e-100) {
    322                     value = 1.0 / value;
    323                     largest = CoinMax(largest, fabs(value));
    324                     sparseFactor_[numberElements] = -value;
    325                     choleskyRow_[numberElements++] = iColumn;
    326                     CoinBigIndex start = columnStart[iColumn];
    327                     CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    328                     for (CoinBigIndex j = start; j < end; j++) {
    329                          choleskyRow_[numberElements] = row[j] + numberTotal;
    330                          sparseFactor_[numberElements++] = element[j];
    331                          largest = CoinMax(largest, fabs(element[j]));
    332                     }
    333                } else {
    334                     sparseFactor_[numberElements] = -1.0e100;
    335                     choleskyRow_[numberElements++] = iColumn;
    336                }
     293    //if (model_->model()->logLevel()&4)
     294    std::cout << "large perturbation " << perturbation << std::endl;
     295#endif
     296    perturbation = sqrt(perturbation);
     297    ;
     298    perturbation = 1.0;
     299  }
     300  // need to recreate every time
     301  int iRow, iColumn;
     302  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
     303  const int *columnLength = model_->clpMatrix()->getVectorLengths();
     304  const int *row = model_->clpMatrix()->getIndices();
     305  const double *element = model_->clpMatrix()->getElements();
     306
     307  CoinBigIndex numberElements = 0;
     308  CoinPackedMatrix *quadratic = NULL;
     309  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(model_->objectiveAsObject()));
     310  if (quadraticObj)
     311    quadratic = quadraticObj->quadraticObjective();
     312  // matrix
     313  if (!quadratic) {
     314    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     315      choleskyStart_[iColumn] = numberElements;
     316      double value = diagonal[iColumn];
     317      if (fabs(value) > 1.0e-100) {
     318        value = 1.0 / value;
     319        largest = CoinMax(largest, fabs(value));
     320        sparseFactor_[numberElements] = -value;
     321        choleskyRow_[numberElements++] = iColumn;
     322        CoinBigIndex start = columnStart[iColumn];
     323        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     324        for (CoinBigIndex j = start; j < end; j++) {
     325          choleskyRow_[numberElements] = row[j] + numberTotal;
     326          sparseFactor_[numberElements++] = element[j];
     327          largest = CoinMax(largest, fabs(element[j]));
     328        }
     329      } else {
     330        sparseFactor_[numberElements] = -1.0e100;
     331        choleskyRow_[numberElements++] = iColumn;
     332      }
     333    }
     334  } else {
     335    // Quadratic
     336    const int *columnQuadratic = quadratic->getIndices();
     337    const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
     338    const int *columnQuadraticLength = quadratic->getVectorLengths();
     339    const double *quadraticElement = quadratic->getElements();
     340    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     341      choleskyStart_[iColumn] = numberElements;
     342      CoinBigIndex savePosition = numberElements;
     343      choleskyRow_[numberElements++] = iColumn;
     344      double value = diagonal[iColumn];
     345      if (fabs(value) > 1.0e-100) {
     346        value = 1.0 / value;
     347        for (CoinBigIndex j = columnQuadraticStart[iColumn];
     348             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     349          int jColumn = columnQuadratic[j];
     350          if (jColumn > iColumn) {
     351            sparseFactor_[numberElements] = -quadraticElement[j];
     352            choleskyRow_[numberElements++] = jColumn;
     353          } else if (iColumn == jColumn) {
     354            value += quadraticElement[j];
    337355          }
    338      } else {
    339           // Quadratic
    340           const int * columnQuadratic = quadratic->getIndices();
    341           const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
    342           const int * columnQuadraticLength = quadratic->getVectorLengths();
    343           const double * quadraticElement = quadratic->getElements();
    344           for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    345                choleskyStart_[iColumn] = numberElements;
    346                CoinBigIndex savePosition = numberElements;
    347                choleskyRow_[numberElements++] = iColumn;
    348                double value = diagonal[iColumn];
    349                if (fabs(value) > 1.0e-100) {
    350                     value = 1.0 / value;
    351                     for (CoinBigIndex j = columnQuadraticStart[iColumn];
    352                               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    353                          int jColumn = columnQuadratic[j];
    354                          if (jColumn > iColumn) {
    355                               sparseFactor_[numberElements] = -quadraticElement[j];
    356                               choleskyRow_[numberElements++] = jColumn;
    357                          } else if (iColumn == jColumn) {
    358                               value += quadraticElement[j];
    359                          }
    360                     }
    361                     largest = CoinMax(largest, fabs(value));
    362                     sparseFactor_[savePosition] = -value;
    363                     CoinBigIndex start = columnStart[iColumn];
    364                     CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    365                     for (CoinBigIndex j = start; j < end; j++) {
    366                          choleskyRow_[numberElements] = row[j] + numberTotal;
    367                          sparseFactor_[numberElements++] = element[j];
    368                          largest = CoinMax(largest, fabs(element[j]));
    369                     }
    370                } else {
    371                     value = 1.0e100;
    372                     sparseFactor_[savePosition] = -value;
    373                }
    374           }
    375      }
    376      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    377           assert (sparseFactor_[choleskyStart_[iColumn]] < 0.0);
    378      }
    379      // slacks
    380      for (iColumn = numberColumns; iColumn < numberTotal; iColumn++) {
    381           choleskyStart_[iColumn] = numberElements;
    382           double value = diagonal[iColumn];
    383           if (fabs(value) > 1.0e-100) {
    384                value = 1.0 / value;
    385                largest = CoinMax(largest, fabs(value));
    386           } else {
    387                value = 1.0e100;
    388           }
    389           sparseFactor_[numberElements] = -value;
    390           choleskyRow_[numberElements++] = iColumn;
    391           choleskyRow_[numberElements] = iColumn - numberColumns + numberTotal;
    392           sparseFactor_[numberElements++] = -1.0;
    393      }
    394      // Finish diagonal
    395      double delta2 = model_->delta(); // add delta*delta to bottom
    396      delta2 *= delta2;
    397      for (iRow = 0; iRow < numberRowsModel; iRow++) {
    398           choleskyStart_[iRow+numberTotal] = numberElements;
    399           choleskyRow_[numberElements] = iRow + numberTotal;
    400           sparseFactor_[numberElements++] = delta2;
    401      }
    402      choleskyStart_[numberRows_] = numberElements;
    403      int i1 = 1;
    404      int i0 = 0;
    405      integerParameters_[1] = 3;
    406      integerParameters_[2] = 3;
    407      integerParameters_[10] = 2;
    408      //integerParameters_[11]=1;
    409      integerParameters_[12] = 2;
    410      // LDLT
    411      integerParameters_[30] = 1;
    412      doubleParameters_[20] = 1.0e100;
    413      double largest2 = largest * 1.0e-20;
    414      largest = CoinMin(largest2, 1.0e-11);
    415      doubleParameters_[10] = CoinMax(1.0e-20, largest);
    416      if (doubleParameters_[10] > 1.0e-3)
    417           integerParameters_[9] = 1;
    418      else
    419           integerParameters_[9] = 0;
     356        }
     357        largest = CoinMax(largest, fabs(value));
     358        sparseFactor_[savePosition] = -value;
     359        CoinBigIndex start = columnStart[iColumn];
     360        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     361        for (CoinBigIndex j = start; j < end; j++) {
     362          choleskyRow_[numberElements] = row[j] + numberTotal;
     363          sparseFactor_[numberElements++] = element[j];
     364          largest = CoinMax(largest, fabs(element[j]));
     365        }
     366      } else {
     367        value = 1.0e100;
     368        sparseFactor_[savePosition] = -value;
     369      }
     370    }
     371  }
     372  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     373    assert(sparseFactor_[choleskyStart_[iColumn]] < 0.0);
     374  }
     375  // slacks
     376  for (iColumn = numberColumns; iColumn < numberTotal; iColumn++) {
     377    choleskyStart_[iColumn] = numberElements;
     378    double value = diagonal[iColumn];
     379    if (fabs(value) > 1.0e-100) {
     380      value = 1.0 / value;
     381      largest = CoinMax(largest, fabs(value));
     382    } else {
     383      value = 1.0e100;
     384    }
     385    sparseFactor_[numberElements] = -value;
     386    choleskyRow_[numberElements++] = iColumn;
     387    choleskyRow_[numberElements] = iColumn - numberColumns + numberTotal;
     388    sparseFactor_[numberElements++] = -1.0;
     389  }
     390  // Finish diagonal
     391  double delta2 = model_->delta(); // add delta*delta to bottom
     392  delta2 *= delta2;
     393  for (iRow = 0; iRow < numberRowsModel; iRow++) {
     394    choleskyStart_[iRow + numberTotal] = numberElements;
     395    choleskyRow_[numberElements] = iRow + numberTotal;
     396    sparseFactor_[numberElements++] = delta2;
     397  }
     398  choleskyStart_[numberRows_] = numberElements;
     399  int i1 = 1;
     400  int i0 = 0;
     401  integerParameters_[1] = 3;
     402  integerParameters_[2] = 3;
     403  integerParameters_[10] = 2;
     404  //integerParameters_[11]=1;
     405  integerParameters_[12] = 2;
     406  // LDLT
     407  integerParameters_[30] = 1;
     408  doubleParameters_[20] = 1.0e100;
     409  double largest2 = largest * 1.0e-20;
     410  largest = CoinMin(largest2, 1.0e-11);
     411  doubleParameters_[10] = CoinMax(1.0e-20, largest);
     412  if (doubleParameters_[10] > 1.0e-3)
     413    integerParameters_[9] = 1;
     414  else
     415    integerParameters_[9] = 0;
    420416#ifndef WSMP
    421      // Set up LDL cutoff
    422      integerParameters_[34] = numberTotal;
    423      doubleParameters_[20] = 1.0e-15;
    424      doubleParameters_[34] = 1.0e-12;
    425      //printf("tol is %g\n",doubleParameters_[10]);
    426      //doubleParameters_[10]=1.0e-17;
    427 #endif
    428      int * rowsDropped2 = new int[numberRows_];
    429      CoinZeroN(rowsDropped2, numberRows_);
    430      F77_FUNC(wssmp,WSSMP)(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
    431            NULL, permute_, permuteInverse_, NULL, &numberRows_, &i1,
    432            NULL, &i0, rowsDropped2, integerParameters_, doubleParameters_);
    433      //std::cout<<"factorization took "<<doubleParameters_[0]<<std::endl;
    434      if (integerParameters_[9]) {
    435           std::cout << "scaling applied" << std::endl;
    436      }
    437      newDropped = integerParameters_[20];
     417  // Set up LDL cutoff
     418  integerParameters_[34] = numberTotal;
     419  doubleParameters_[20] = 1.0e-15;
     420  doubleParameters_[34] = 1.0e-12;
     421  //printf("tol is %g\n",doubleParameters_[10]);
     422  //doubleParameters_[10]=1.0e-17;
     423#endif
     424  int *rowsDropped2 = new int[numberRows_];
     425  CoinZeroN(rowsDropped2, numberRows_);
     426  F77_FUNC(wssmp, WSSMP)
     427  (&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
     428    NULL, permute_, permuteInverse_, NULL, &numberRows_, &i1,
     429    NULL, &i0, rowsDropped2, integerParameters_, doubleParameters_);
     430  //std::cout<<"factorization took "<<doubleParameters_[0]<<std::endl;
     431  if (integerParameters_[9]) {
     432    std::cout << "scaling applied" << std::endl;
     433  }
     434  newDropped = integerParameters_[20];
    438435#if 1
    439      // Should save adjustments in ..R_
    440      int n1 = 0, n2 = 0;
    441      double * primalR = model_->primalR();
    442      double * dualR = model_->dualR();
    443      for (iRow = 0; iRow < numberTotal; iRow++) {
    444           if (rowsDropped2[iRow]) {
    445                n1++;
    446                //printf("row region1 %d dropped\n",iRow);
    447                //rowsDropped_[iRow]=1;
    448                rowsDropped_[iRow] = 0;
    449                primalR[iRow] = doubleParameters_[20];
    450           } else {
    451                rowsDropped_[iRow] = 0;
    452                primalR[iRow] = 0.0;
    453           }
    454      }
    455      for (; iRow < numberRows_; iRow++) {
    456           if (rowsDropped2[iRow]) {
    457                n2++;
    458                //printf("row region2 %d dropped\n",iRow);
    459                //rowsDropped_[iRow]=1;
    460                rowsDropped_[iRow] = 0;
    461                dualR[iRow-numberTotal] = doubleParameters_[34];
    462           } else {
    463                rowsDropped_[iRow] = 0;
    464                dualR[iRow-numberTotal] = 0.0;
    465           }
    466      }
    467      //printf("%d rows dropped in region1, %d in region2\n",n1,n2);
    468 #endif
    469      delete [] rowsDropped2;
    470      //if (integerParameters_[20])
    471      //std::cout<<integerParameters_[20]<<" rows dropped"<<std::endl;
    472      largest = doubleParameters_[3];
    473      smallest = doubleParameters_[4];
    474      if (model_->messageHandler()->logLevel() > 1)
    475           std::cout << "Cholesky - largest " << largest << " smallest " << smallest << std::endl;
    476      choleskyCondition_ = largest / smallest;
    477      if (integerParameters_[63] < 0)
    478           return -1; // out of memory
    479      status_ = 0;
    480      return 0;
     436  // Should save adjustments in ..R_
     437  int n1 = 0, n2 = 0;
     438  double *primalR = model_->primalR();
     439  double *dualR = model_->dualR();
     440  for (iRow = 0; iRow < numberTotal; iRow++) {
     441    if (rowsDropped2[iRow]) {
     442      n1++;
     443      //printf("row region1 %d dropped\n",iRow);
     444      //rowsDropped_[iRow]=1;
     445      rowsDropped_[iRow] = 0;
     446      primalR[iRow] = doubleParameters_[20];
     447    } else {
     448      rowsDropped_[iRow] = 0;
     449      primalR[iRow] = 0.0;
     450    }
     451  }
     452  for (; iRow < numberRows_; iRow++) {
     453    if (rowsDropped2[iRow]) {
     454      n2++;
     455      //printf("row region2 %d dropped\n",iRow);
     456      //rowsDropped_[iRow]=1;
     457      rowsDropped_[iRow] = 0;
     458      dualR[iRow - numberTotal] = doubleParameters_[34];
     459    } else {
     460      rowsDropped_[iRow] = 0;
     461      dualR[iRow - numberTotal] = 0.0;
     462    }
     463  }
     464  //printf("%d rows dropped in region1, %d in region2\n",n1,n2);
     465#endif
     466  delete[] rowsDropped2;
     467  //if (integerParameters_[20])
     468  //std::cout<<integerParameters_[20]<<" rows dropped"<<std::endl;
     469  largest = doubleParameters_[3];
     470  smallest = doubleParameters_[4];
     471  if (model_->messageHandler()->logLevel() > 1)
     472    std::cout << "Cholesky - largest " << largest << " smallest " << smallest << std::endl;
     473  choleskyCondition_ = largest / smallest;
     474  if (integerParameters_[63] < 0)
     475    return -1; // out of memory
     476  status_ = 0;
     477  return 0;
    481478}
    482479/* Uses factorization to solve. */
    483 void
    484 ClpCholeskyWssmpKKT::solve (double * region)
    485 {
    486      abort();
     480void ClpCholeskyWssmpKKT::solve(double *region)
     481{
     482  abort();
    487483}
    488484/* Uses factorization to solve. */
    489 void
    490 ClpCholeskyWssmpKKT::solveKKT (double * region1, double * region2, const double * diagonal,
    491                                double diagonalScaleFactor)
    492 {
    493      int numberRowsModel = model_->numberRows();
    494      int numberColumns = model_->numberColumns();
    495      int numberTotal = numberColumns + numberRowsModel;
    496      double * array = new double [numberRows_];
    497      CoinMemcpyN(region1, numberTotal, array);
    498      CoinMemcpyN(region2, numberRowsModel, array + numberTotal);
    499      int i1 = 1;
    500      int i0 = 0;
    501      integerParameters_[1] = 4;
    502      integerParameters_[2] = 4;
     485void ClpCholeskyWssmpKKT::solveKKT(double *region1, double *region2, const double *diagonal,
     486  double diagonalScaleFactor)
     487{
     488  int numberRowsModel = model_->numberRows();
     489  int numberColumns = model_->numberColumns();
     490  int numberTotal = numberColumns + numberRowsModel;
     491  double *array = new double[numberRows_];
     492  CoinMemcpyN(region1, numberTotal, array);
     493  CoinMemcpyN(region2, numberRowsModel, array + numberTotal);
     494  int i1 = 1;
     495  int i0 = 0;
     496  integerParameters_[1] = 4;
     497  integerParameters_[2] = 4;
    503498#if 0
    504499     integerParameters_[5] = 3;
     
    506501     integerParameters_[6] = 6;
    507502#endif
    508      F77_FUNC(wssmp,WSSMP)(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
    509            NULL, permute_, permuteInverse_, array, &numberRows_, &i1,
    510            NULL, &i0, NULL, integerParameters_, doubleParameters_);
     503  F77_FUNC(wssmp, WSSMP)
     504  (&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
     505    NULL, permute_, permuteInverse_, array, &numberRows_, &i1,
     506    NULL, &i0, NULL, integerParameters_, doubleParameters_);
    511507#if 0
    512508     int iRow;
     
    522518     }
    523519#endif
    524      CoinMemcpyN(array + numberTotal, numberRowsModel, region2);
     520  CoinMemcpyN(array + numberTotal, numberRowsModel, region2);
    525521#if 1
    526      CoinMemcpyN(array, numberTotal, region1);
     522  CoinMemcpyN(array, numberTotal, region1);
    527523#else
    528      multiplyAdd(region2, numberRowsModel, -1.0, array + numberColumns, 0.0);
    529      CoinZeroN(array, numberColumns);
    530      model_->clpMatrix()->transposeTimes(1.0, region2, array);
    531      for (int iColumn = 0; iColumn < numberTotal; iColumn++)
    532           region1[iColumn] = diagonal[iColumn] * (array[iColumn] - region1[iColumn]);
    533 #endif
    534      delete [] array;
     524  multiplyAdd(region2, numberRowsModel, -1.0, array + numberColumns, 0.0);
     525  CoinZeroN(array, numberColumns);
     526  model_->clpMatrix()->transposeTimes(1.0, region2, array);
     527  for (int iColumn = 0; iColumn < numberTotal; iColumn++)
     528    region1[iColumn] = diagonal[iColumn] * (array[iColumn] - region1[iColumn]);
     529#endif
     530  delete[] array;
    535531#if 0
    536532     if (integerParameters_[5]) {
     
    540536#endif
    541537}
     538
     539/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     540*/
Note: See TracChangeset for help on using the changeset viewer.