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/ClpCholeskyUfl.cpp

    r1723 r2385  
    3131// Default Constructor
    3232//-------------------------------------------------------------------
    33 ClpCholeskyUfl::ClpCholeskyUfl (int denseThreshold)
    34      : ClpCholeskyBase(denseThreshold)
    35 {
    36      type_ = 14;
    37      L_ = NULL;
    38      c_ = NULL;
    39      
     33ClpCholeskyUfl::ClpCholeskyUfl(int denseThreshold)
     34  : ClpCholeskyBase(denseThreshold)
     35{
     36  type_ = 14;
     37  L_ = NULL;
     38  c_ = NULL;
     39
    4040#ifdef COIN_HAS_CHOLMOD
    41      c_ = (cholmod_common*) malloc(sizeof(cholmod_common));
    42      cholmod_start (c_) ;
    43      // Can't use supernodal as may not be positive definite
    44      c_->supernodal = 0;
     41  c_ = (cholmod_common *)malloc(sizeof(cholmod_common));
     42  cholmod_start(c_);
     43  // Can't use supernodal as may not be positive definite
     44  c_->supernodal = 0;
    4545#endif
    4646}
     
    4949// Copy constructor
    5050//-------------------------------------------------------------------
    51 ClpCholeskyUfl::ClpCholeskyUfl (const ClpCholeskyUfl & rhs)
    52      : ClpCholeskyBase(rhs)
    53 {
    54      abort();
    55 }
    56 
     51ClpCholeskyUfl::ClpCholeskyUfl(const ClpCholeskyUfl &rhs)
     52  : ClpCholeskyBase(rhs)
     53{
     54  abort();
     55}
    5756
    5857//-------------------------------------------------------------------
    5958// Destructor
    6059//-------------------------------------------------------------------
    61 ClpCholeskyUfl::~ClpCholeskyUfl ()
     60ClpCholeskyUfl::~ClpCholeskyUfl()
    6261{
    6362#ifdef COIN_HAS_CHOLMOD
    64      cholmod_free_factor (&L_, c_) ;
    65      cholmod_finish (c_) ;
    66      free(c_);
     63  cholmod_free_factor(&L_, c_);
     64  cholmod_finish(c_);
     65  free(c_);
    6766#endif
    6867}
     
    7271//-------------------------------------------------------------------
    7372ClpCholeskyUfl &
    74 ClpCholeskyUfl::operator=(const ClpCholeskyUfl& rhs)
    75 {
    76      if (this != &rhs) {
    77           ClpCholeskyBase::operator=(rhs);
    78           abort();
    79      }
    80      return *this;
     73ClpCholeskyUfl::operator=(const ClpCholeskyUfl &rhs)
     74{
     75  if (this != &rhs) {
     76    ClpCholeskyBase::operator=(rhs);
     77    abort();
     78  }
     79  return *this;
    8180}
    8281
     
    8483// Clone
    8584//-------------------------------------------------------------------
    86 ClpCholeskyBase * ClpCholeskyUfl::clone() const
    87 {
    88      return new ClpCholeskyUfl(*this);
     85ClpCholeskyBase *ClpCholeskyUfl::clone() const
     86{
     87  return new ClpCholeskyUfl(*this);
    8988}
    9089
    9190#ifndef COIN_HAS_CHOLMOD
    9291/* Orders rows and saves pointer to matrix.and model */
    93 int
    94 ClpCholeskyUfl::order(ClpInterior * model)
    95 {
    96      int iRow;
    97      model_ = model;
    98      if (preOrder(false, true, doKKT_))
    99           return -1;
    100      permuteInverse_ = new int [numberRows_];
    101      permute_ = new int[numberRows_];
    102      double Control[AMD_CONTROL];
    103      double Info[AMD_INFO];
    104 
    105      amd_defaults(Control);
    106      //amd_control(Control);
    107 
    108      int returnCode = amd_order(numberRows_, choleskyStart_, choleskyRow_,
    109                                 permute_, Control, Info);
    110      delete [] choleskyRow_;
    111      choleskyRow_ = NULL;
    112      delete [] choleskyStart_;
    113      choleskyStart_ = NULL;
    114      //amd_info(Info);
    115 
    116      if (returnCode != AMD_OK) {
    117           std::cout << "AMD ordering failed" << std::endl;
    118           return 1;
    119      }
    120      for (iRow = 0; iRow < numberRows_; iRow++) {
    121           permuteInverse_[permute_[iRow]] = iRow;
    122      }
    123      return 0;
     92int ClpCholeskyUfl::order(ClpInterior *model)
     93{
     94  int iRow;
     95  model_ = model;
     96  if (preOrder(false, true, doKKT_))
     97    return -1;
     98  permuteInverse_ = new int[numberRows_];
     99  permute_ = new int[numberRows_];
     100  double Control[AMD_CONTROL];
     101  double Info[AMD_INFO];
     102
     103  amd_defaults(Control);
     104  //amd_control(Control);
     105
     106  int returnCode = amd_order(numberRows_, choleskyStart_, choleskyRow_,
     107    permute_, Control, Info);
     108  delete[] choleskyRow_;
     109  choleskyRow_ = NULL;
     110  delete[] choleskyStart_;
     111  choleskyStart_ = NULL;
     112  //amd_info(Info);
     113
     114  if (returnCode != AMD_OK) {
     115    std::cout << "AMD ordering failed" << std::endl;
     116    return 1;
     117  }
     118  for (iRow = 0; iRow < numberRows_; iRow++) {
     119    permuteInverse_[permute_[iRow]] = iRow;
     120  }
     121  return 0;
    124122}
    125123#else
    126124/* Orders rows and saves pointer to matrix.and model */
    127 int
    128 ClpCholeskyUfl::order(ClpInterior * model)
    129 {
    130      numberRows_ = model->numberRows();
    131      if (doKKT_) {
    132           numberRows_ += numberRows_ + model->numberColumns();
    133           printf("finish coding UFL KKT!\n");
    134           abort();
    135      }
    136      rowsDropped_ = new char [numberRows_];
    137      memset(rowsDropped_, 0, numberRows_);
    138      numberRowsDropped_ = 0;
    139      model_ = model;
    140      rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
    141      // Space for starts
    142      choleskyStart_ = new CoinBigIndex[numberRows_+1];
    143      const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    144      const int * columnLength = model_->clpMatrix()->getVectorLengths();
    145      const int * row = model_->clpMatrix()->getIndices();
    146      const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
    147      const int * rowLength = rowCopy_->getVectorLengths();
    148      const int * column = rowCopy_->getIndices();
    149      // We need two arrays for counts
    150      int * which = new int [numberRows_];
    151      int * used = new int[numberRows_+1];
    152      CoinZeroN(used, numberRows_);
    153      int iRow;
    154      sizeFactor_ = 0;
    155      for (iRow = 0; iRow < numberRows_; iRow++) {
    156           int number = 1;
    157           // make sure diagonal exists
    158           which[0] = iRow;
    159           used[iRow] = 1;
    160           if (!rowsDropped_[iRow]) {
    161                CoinBigIndex startRow = rowStart[iRow];
    162                CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
    163                for (CoinBigIndex k = startRow; k < endRow; k++) {
    164                     int iColumn = column[k];
    165                     CoinBigIndex start = columnStart[iColumn];
    166                     CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    167                     for (CoinBigIndex j = start; j < end; j++) {
    168                          int jRow = row[j];
    169                          if (jRow >= iRow && !rowsDropped_[jRow]) {
    170                               if (!used[jRow]) {
    171                                    used[jRow] = 1;
    172                                    which[number++] = jRow;
    173                               }
    174                          }
    175                     }
    176                }
    177                sizeFactor_ += number;
    178                int j;
    179                for (j = 0; j < number; j++)
    180                     used[which[j]] = 0;
     125int ClpCholeskyUfl::order(ClpInterior *model)
     126{
     127  numberRows_ = model->numberRows();
     128  if (doKKT_) {
     129    numberRows_ += numberRows_ + model->numberColumns();
     130    printf("finish coding UFL KKT!\n");
     131    abort();
     132  }
     133  rowsDropped_ = new char[numberRows_];
     134  memset(rowsDropped_, 0, numberRows_);
     135  numberRowsDropped_ = 0;
     136  model_ = model;
     137  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
     138  // Space for starts
     139  choleskyStart_ = new CoinBigIndex[numberRows_ + 1];
     140  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
     141  const int *columnLength = model_->clpMatrix()->getVectorLengths();
     142  const int *row = model_->clpMatrix()->getIndices();
     143  const CoinBigIndex *rowStart = rowCopy_->getVectorStarts();
     144  const int *rowLength = rowCopy_->getVectorLengths();
     145  const int *column = rowCopy_->getIndices();
     146  // We need two arrays for counts
     147  int *which = new int[numberRows_];
     148  int *used = new int[numberRows_ + 1];
     149  CoinZeroN(used, numberRows_);
     150  int iRow;
     151  sizeFactor_ = 0;
     152  for (iRow = 0; iRow < numberRows_; iRow++) {
     153    int number = 1;
     154    // make sure diagonal exists
     155    which[0] = iRow;
     156    used[iRow] = 1;
     157    if (!rowsDropped_[iRow]) {
     158      CoinBigIndex startRow = rowStart[iRow];
     159      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
     160      for (CoinBigIndex k = startRow; k < endRow; k++) {
     161        int iColumn = column[k];
     162        CoinBigIndex start = columnStart[iColumn];
     163        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     164        for (CoinBigIndex j = start; j < end; j++) {
     165          int jRow = row[j];
     166          if (jRow >= iRow && !rowsDropped_[jRow]) {
     167            if (!used[jRow]) {
     168              used[jRow] = 1;
     169              which[number++] = jRow;
     170            }
    181171          }
    182      }
    183      delete [] which;
    184      // Now we have size - create arrays and fill in
    185      try {
    186           choleskyRow_ = new int [sizeFactor_];
    187      } catch (...) {
    188           // no memory
    189           delete [] choleskyStart_;
    190           choleskyStart_ = NULL;
    191           return -1;
    192      }
    193      try {
    194           sparseFactor_ = new double[sizeFactor_];
    195      } catch (...) {
    196           // no memory
    197           delete [] choleskyRow_;
    198           choleskyRow_ = NULL;
    199           delete [] choleskyStart_;
    200           choleskyStart_ = NULL;
    201           return -1;
    202      }
    203 
    204      sizeFactor_ = 0;
    205      which = choleskyRow_;
    206      for (iRow = 0; iRow < numberRows_; iRow++) {
    207           int number = 1;
    208           // make sure diagonal exists
    209           which[0] = iRow;
    210           used[iRow] = 1;
    211           choleskyStart_[iRow] = sizeFactor_;
    212           if (!rowsDropped_[iRow]) {
    213                CoinBigIndex startRow = rowStart[iRow];
    214                CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
    215                for (CoinBigIndex k = startRow; k < endRow; k++) {
    216                     int iColumn = column[k];
    217                     CoinBigIndex start = columnStart[iColumn];
    218                     CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    219                     for (CoinBigIndex j = start; j < end; j++) {
    220                          int jRow = row[j];
    221                          if (jRow >= iRow && !rowsDropped_[jRow]) {
    222                               if (!used[jRow]) {
    223                                    used[jRow] = 1;
    224                                    which[number++] = jRow;
    225                               }
    226                          }
    227                     }
    228                }
    229                sizeFactor_ += number;
    230                int j;
    231                for (j = 0; j < number; j++)
    232                     used[which[j]] = 0;
    233                // Sort
    234                std::sort(which, which + number);
    235                // move which on
    236                which += number;
     172        }
     173      }
     174      sizeFactor_ += number;
     175      int j;
     176      for (j = 0; j < number; j++)
     177        used[which[j]] = 0;
     178    }
     179  }
     180  delete[] which;
     181  // Now we have size - create arrays and fill in
     182  try {
     183    choleskyRow_ = new int[sizeFactor_];
     184  } catch (...) {
     185    // no memory
     186    delete[] choleskyStart_;
     187    choleskyStart_ = NULL;
     188    return -1;
     189  }
     190  try {
     191    sparseFactor_ = new double[sizeFactor_];
     192  } catch (...) {
     193    // no memory
     194    delete[] choleskyRow_;
     195    choleskyRow_ = NULL;
     196    delete[] choleskyStart_;
     197    choleskyStart_ = NULL;
     198    return -1;
     199  }
     200
     201  sizeFactor_ = 0;
     202  which = choleskyRow_;
     203  for (iRow = 0; iRow < numberRows_; iRow++) {
     204    int number = 1;
     205    // make sure diagonal exists
     206    which[0] = iRow;
     207    used[iRow] = 1;
     208    choleskyStart_[iRow] = sizeFactor_;
     209    if (!rowsDropped_[iRow]) {
     210      CoinBigIndex startRow = rowStart[iRow];
     211      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
     212      for (CoinBigIndex k = startRow; k < endRow; k++) {
     213        int iColumn = column[k];
     214        CoinBigIndex start = columnStart[iColumn];
     215        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     216        for (CoinBigIndex j = start; j < end; j++) {
     217          int jRow = row[j];
     218          if (jRow >= iRow && !rowsDropped_[jRow]) {
     219            if (!used[jRow]) {
     220              used[jRow] = 1;
     221              which[number++] = jRow;
     222            }
    237223          }
    238      }
    239      choleskyStart_[numberRows_] = sizeFactor_;
    240      delete [] used;
    241      permuteInverse_ = new int [numberRows_];
    242      permute_ = new int[numberRows_];
    243      cholmod_sparse A ;
    244      A.nrow = numberRows_;
    245      A.ncol = numberRows_;
    246      A.nzmax = choleskyStart_[numberRows_];
    247      A.p = choleskyStart_;
    248      A.i = choleskyRow_;
    249      A.x = NULL;
    250      A.stype = -1;
    251      A.itype = CHOLMOD_INT;
    252      A.xtype = CHOLMOD_PATTERN;
    253      A.dtype = CHOLMOD_DOUBLE;
    254      A.sorted = 1;
    255      A.packed = 1;
    256      c_->nmethods = 9;
    257      c_->postorder = true;
    258      //c_->dbound=1.0e-20;
    259      L_ = cholmod_analyze (&A, c_) ;
    260      if (c_->status) {
    261        COIN_DETAIL_PRINT(std::cout << "CHOLMOD ordering failed" << std::endl);
    262           return 1;
    263      } else {
    264        COIN_DETAIL_PRINT(printf("%g nonzeros, flop count %g\n", c_->lnz, c_->fl));
    265      }
    266      for (iRow = 0; iRow < numberRows_; iRow++) {
    267           permuteInverse_[iRow] = iRow;
    268           permute_[iRow] = iRow;
    269      }
    270      return 0;
     224        }
     225      }
     226      sizeFactor_ += number;
     227      int j;
     228      for (j = 0; j < number; j++)
     229        used[which[j]] = 0;
     230      // Sort
     231      std::sort(which, which + number);
     232      // move which on
     233      which += number;
     234    }
     235  }
     236  choleskyStart_[numberRows_] = sizeFactor_;
     237  delete[] used;
     238  permuteInverse_ = new int[numberRows_];
     239  permute_ = new int[numberRows_];
     240  cholmod_sparse A;
     241  A.nrow = numberRows_;
     242  A.ncol = numberRows_;
     243  A.nzmax = choleskyStart_[numberRows_];
     244  A.p = choleskyStart_;
     245  A.i = choleskyRow_;
     246  A.x = NULL;
     247  A.stype = -1;
     248  A.itype = CHOLMOD_INT;
     249  A.xtype = CHOLMOD_PATTERN;
     250  A.dtype = CHOLMOD_DOUBLE;
     251  A.sorted = 1;
     252  A.packed = 1;
     253  c_->nmethods = 9;
     254  c_->postorder = true;
     255  //c_->dbound=1.0e-20;
     256  L_ = cholmod_analyze(&A, c_);
     257  if (c_->status) {
     258    COIN_DETAIL_PRINT(std::cout << "CHOLMOD ordering failed" << std::endl);
     259    return 1;
     260  } else {
     261    COIN_DETAIL_PRINT(printf("%g nonzeros, flop count %g\n", c_->lnz, c_->fl));
     262  }
     263  for (iRow = 0; iRow < numberRows_; iRow++) {
     264    permuteInverse_[iRow] = iRow;
     265    permute_[iRow] = iRow;
     266  }
     267  return 0;
    271268}
    272269#endif
     
    276273   user must provide factorize and solve.  Otherwise the default factorization is used
    277274   returns non-zero if not enough memory */
    278 int
    279 ClpCholeskyUfl::symbolic()
     275int ClpCholeskyUfl::symbolic()
    280276{
    281277#ifdef COIN_HAS_CHOLMOD
    282      return 0;
    283 #else
    284      return ClpCholeskyBase::symbolic();
     278  return 0;
     279#else
     280  return ClpCholeskyBase::symbolic();
    285281#endif
    286282}
     
    288284#ifdef COIN_HAS_CHOLMOD
    289285/* Factorize - filling in rowsDropped and returning number dropped */
    290 int
    291 ClpCholeskyUfl::factorize(const double * diagonal, int * rowsDropped)
    292 {
    293      const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    294      const int * columnLength = model_->clpMatrix()->getVectorLengths();
    295      const int * row = model_->clpMatrix()->getIndices();
    296      const double * element = model_->clpMatrix()->getElements();
    297      const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
    298      const int * rowLength = rowCopy_->getVectorLengths();
    299      const int * column = rowCopy_->getIndices();
    300      const double * elementByRow = rowCopy_->getElements();
    301      int numberColumns = model_->clpMatrix()->getNumCols();
    302      int iRow;
    303      double * work = new double[numberRows_];
    304      CoinZeroN(work, numberRows_);
    305      const double * diagonalSlack = diagonal + numberColumns;
    306      int newDropped = 0;
    307      double largest;
    308      //double smallest;
    309      //perturbation
    310      double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
    311      perturbation = 0.0;
    312      perturbation = perturbation * perturbation;
    313      if (perturbation > 1.0) {
     286int ClpCholeskyUfl::factorize(const double *diagonal, int *rowsDropped)
     287{
     288  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
     289  const int *columnLength = model_->clpMatrix()->getVectorLengths();
     290  const int *row = model_->clpMatrix()->getIndices();
     291  const double *element = model_->clpMatrix()->getElements();
     292  const CoinBigIndex *rowStart = rowCopy_->getVectorStarts();
     293  const int *rowLength = rowCopy_->getVectorLengths();
     294  const int *column = rowCopy_->getIndices();
     295  const double *elementByRow = rowCopy_->getElements();
     296  int numberColumns = model_->clpMatrix()->getNumCols();
     297  int iRow;
     298  double *work = new double[numberRows_];
     299  CoinZeroN(work, numberRows_);
     300  const double *diagonalSlack = diagonal + numberColumns;
     301  int newDropped = 0;
     302  double largest;
     303  //double smallest;
     304  //perturbation
     305  double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
     306  perturbation = 0.0;
     307  perturbation = perturbation * perturbation;
     308  if (perturbation > 1.0) {
    314309#ifdef COIN_DEVELOP
    315           //if (model_->model()->logLevel()&4)
    316           std::cout << "large perturbation " << perturbation << std::endl;
    317 #endif
    318           perturbation = sqrt(perturbation);;
    319           perturbation = 1.0;
    320      }
    321      double delta2 = model_->delta(); // add delta*delta to diagonal
    322      delta2 *= delta2;
    323      for (iRow = 0; iRow < numberRows_; iRow++) {
    324           double * put = sparseFactor_ + choleskyStart_[iRow];
    325           int * which = choleskyRow_ + choleskyStart_[iRow];
    326           int number = choleskyStart_[iRow+1] - choleskyStart_[iRow];
    327           if (!rowLength[iRow])
    328                rowsDropped_[iRow] = 1;
    329           if (!rowsDropped_[iRow]) {
    330                CoinBigIndex startRow = rowStart[iRow];
    331                CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
    332                work[iRow] = diagonalSlack[iRow] + delta2;
    333                for (CoinBigIndex k = startRow; k < endRow; k++) {
    334                     int iColumn = column[k];
    335                     if (!whichDense_ || !whichDense_[iColumn]) {
    336                          CoinBigIndex start = columnStart[iColumn];
    337                          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    338                          double multiplier = diagonal[iColumn] * elementByRow[k];
    339                          for (CoinBigIndex j = start; j < end; j++) {
    340                               int jRow = row[j];
    341                               if (jRow >= iRow && !rowsDropped_[jRow]) {
    342                                    double value = element[j] * multiplier;
    343                                    work[jRow] += value;
    344                               }
    345                          }
    346                     }
    347                }
    348                int j;
    349                for (j = 0; j < number; j++) {
    350                     int jRow = which[j];
    351                     put[j] = work[jRow];
    352                     work[jRow] = 0.0;
    353                }
    354           } else {
    355                // dropped
    356                int j;
    357                for (j = 1; j < number; j++) {
    358                     put[j] = 0.0;
    359                }
    360                put[0] = 1.0;
     310    //if (model_->model()->logLevel()&4)
     311    std::cout << "large perturbation " << perturbation << std::endl;
     312#endif
     313    perturbation = sqrt(perturbation);
     314    ;
     315    perturbation = 1.0;
     316  }
     317  double delta2 = model_->delta(); // add delta*delta to diagonal
     318  delta2 *= delta2;
     319  for (iRow = 0; iRow < numberRows_; iRow++) {
     320    double *put = sparseFactor_ + choleskyStart_[iRow];
     321    int *which = choleskyRow_ + choleskyStart_[iRow];
     322    int number = choleskyStart_[iRow + 1] - choleskyStart_[iRow];
     323    if (!rowLength[iRow])
     324      rowsDropped_[iRow] = 1;
     325    if (!rowsDropped_[iRow]) {
     326      CoinBigIndex startRow = rowStart[iRow];
     327      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
     328      work[iRow] = diagonalSlack[iRow] + delta2;
     329      for (CoinBigIndex k = startRow; k < endRow; k++) {
     330        int iColumn = column[k];
     331        if (!whichDense_ || !whichDense_[iColumn]) {
     332          CoinBigIndex start = columnStart[iColumn];
     333          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     334          double multiplier = diagonal[iColumn] * elementByRow[k];
     335          for (CoinBigIndex j = start; j < end; j++) {
     336            int jRow = row[j];
     337            if (jRow >= iRow && !rowsDropped_[jRow]) {
     338              double value = element[j] * multiplier;
     339              work[jRow] += value;
     340            }
    361341          }
    362      }
    363      //check sizes
    364      double largest2 = maximumAbsElement(sparseFactor_, sizeFactor_);
    365      largest2 *= 1.0e-20;
    366      largest = CoinMin(largest2, 1.0e-11);
    367      int numberDroppedBefore = 0;
    368      for (iRow = 0; iRow < numberRows_; iRow++) {
    369           int dropped = rowsDropped_[iRow];
    370           // Move to int array
    371           rowsDropped[iRow] = dropped;
    372           if (!dropped) {
    373                CoinBigIndex start = choleskyStart_[iRow];
    374                double diagonal = sparseFactor_[start];
    375                if (diagonal > largest2) {
    376                     sparseFactor_[start] = CoinMax(diagonal, 1.0e-10);
    377                } else {
    378                     sparseFactor_[start] = CoinMax(diagonal, 1.0e-10);
    379                     rowsDropped[iRow] = 2;
    380                     numberDroppedBefore++;
    381                }
    382           }
    383      }
    384      delete [] work;
    385      cholmod_sparse A ;
    386      A.nrow = numberRows_;
    387      A.ncol = numberRows_;
    388      A.nzmax = choleskyStart_[numberRows_];
    389      A.p = choleskyStart_;
    390      A.i = choleskyRow_;
    391      A.x = sparseFactor_;
    392      A.stype = -1;
    393      A.itype = CHOLMOD_INT;
    394      A.xtype = CHOLMOD_REAL;
    395      A.dtype = CHOLMOD_DOUBLE;
    396      A.sorted = 1;
    397      A.packed = 1;
    398      cholmod_factorize (&A, L_, c_) ;               /* factorize */
    399      choleskyCondition_ = 1.0;
    400      bool cleanCholesky;
    401      if (model_->numberIterations() < 2000)
    402           cleanCholesky = true;
    403      else
    404           cleanCholesky = false;
    405      if (cleanCholesky) {
    406           //drop fresh makes some formADAT easier
    407           //int oldDropped=numberRowsDropped_;
    408           if (newDropped || numberRowsDropped_) {
    409                //std::cout <<"Rank "<<numberRows_-newDropped<<" ( "<<
    410                //  newDropped<<" dropped)";
    411                //if (newDropped>oldDropped)
    412                //std::cout<<" ( "<<newDropped-oldDropped<<" dropped this time)";
    413                //std::cout<<std::endl;
    414                newDropped = 0;
    415                for (int i = 0; i < numberRows_; i++) {
    416                     int dropped = rowsDropped[i];
    417                     rowsDropped_[i] = (char)dropped;
    418                     if (dropped == 2) {
    419                          //dropped this time
    420                          rowsDropped[newDropped++] = i;
    421                          rowsDropped_[i] = 0;
    422                     }
    423                }
    424                numberRowsDropped_ = newDropped;
    425                newDropped = -(2 + newDropped);
    426           }
    427      } else {
    428           if (newDropped) {
    429                newDropped = 0;
    430                for (int i = 0; i < numberRows_; i++) {
    431                     int dropped = rowsDropped[i];
    432                     rowsDropped_[i] = (char)dropped;
    433                     if (dropped == 2) {
    434                          //dropped this time
    435                          rowsDropped[newDropped++] = i;
    436                          rowsDropped_[i] = 1;
    437                     }
    438                }
    439           }
    440           numberRowsDropped_ += newDropped;
    441           if (numberRowsDropped_ && 0) {
    442                std::cout << "Rank " << numberRows_ - numberRowsDropped_ << " ( " <<
    443                          numberRowsDropped_ << " dropped)";
    444                if (newDropped) {
    445                     std::cout << " ( " << newDropped << " dropped this time)";
    446                }
    447                std::cout << std::endl;
    448           }
    449      }
    450      status_ = 0;
    451      return newDropped;
     342        }
     343      }
     344      int j;
     345      for (j = 0; j < number; j++) {
     346        int jRow = which[j];
     347        put[j] = work[jRow];
     348        work[jRow] = 0.0;
     349      }
     350    } else {
     351      // dropped
     352      int j;
     353      for (j = 1; j < number; j++) {
     354        put[j] = 0.0;
     355      }
     356      put[0] = 1.0;
     357    }
     358  }
     359  //check sizes
     360  double largest2 = maximumAbsElement(sparseFactor_, sizeFactor_);
     361  largest2 *= 1.0e-20;
     362  largest = CoinMin(largest2, 1.0e-11);
     363  int numberDroppedBefore = 0;
     364  for (iRow = 0; iRow < numberRows_; iRow++) {
     365    int dropped = rowsDropped_[iRow];
     366    // Move to int array
     367    rowsDropped[iRow] = dropped;
     368    if (!dropped) {
     369      CoinBigIndex start = choleskyStart_[iRow];
     370      double diagonal = sparseFactor_[start];
     371      if (diagonal > largest2) {
     372        sparseFactor_[start] = CoinMax(diagonal, 1.0e-10);
     373      } else {
     374        sparseFactor_[start] = CoinMax(diagonal, 1.0e-10);
     375        rowsDropped[iRow] = 2;
     376        numberDroppedBefore++;
     377      }
     378    }
     379  }
     380  delete[] work;
     381  cholmod_sparse A;
     382  A.nrow = numberRows_;
     383  A.ncol = numberRows_;
     384  A.nzmax = choleskyStart_[numberRows_];
     385  A.p = choleskyStart_;
     386  A.i = choleskyRow_;
     387  A.x = sparseFactor_;
     388  A.stype = -1;
     389  A.itype = CHOLMOD_INT;
     390  A.xtype = CHOLMOD_REAL;
     391  A.dtype = CHOLMOD_DOUBLE;
     392  A.sorted = 1;
     393  A.packed = 1;
     394  cholmod_factorize(&A, L_, c_); /* factorize */
     395  choleskyCondition_ = 1.0;
     396  bool cleanCholesky;
     397  if (model_->numberIterations() < 2000)
     398    cleanCholesky = true;
     399  else
     400    cleanCholesky = false;
     401  if (cleanCholesky) {
     402    //drop fresh makes some formADAT easier
     403    //int oldDropped=numberRowsDropped_;
     404    if (newDropped || numberRowsDropped_) {
     405      //std::cout <<"Rank "<<numberRows_-newDropped<<" ( "<<
     406      //  newDropped<<" dropped)";
     407      //if (newDropped>oldDropped)
     408      //std::cout<<" ( "<<newDropped-oldDropped<<" dropped this time)";
     409      //std::cout<<std::endl;
     410      newDropped = 0;
     411      for (int i = 0; i < numberRows_; i++) {
     412        int dropped = rowsDropped[i];
     413        rowsDropped_[i] = (char)dropped;
     414        if (dropped == 2) {
     415          //dropped this time
     416          rowsDropped[newDropped++] = i;
     417          rowsDropped_[i] = 0;
     418        }
     419      }
     420      numberRowsDropped_ = newDropped;
     421      newDropped = -(2 + newDropped);
     422    }
     423  } else {
     424    if (newDropped) {
     425      newDropped = 0;
     426      for (int i = 0; i < numberRows_; i++) {
     427        int dropped = rowsDropped[i];
     428        rowsDropped_[i] = (char)dropped;
     429        if (dropped == 2) {
     430          //dropped this time
     431          rowsDropped[newDropped++] = i;
     432          rowsDropped_[i] = 1;
     433        }
     434      }
     435    }
     436    numberRowsDropped_ += newDropped;
     437    if (numberRowsDropped_ && 0) {
     438      std::cout << "Rank " << numberRows_ - numberRowsDropped_ << " ( " << numberRowsDropped_ << " dropped)";
     439      if (newDropped) {
     440        std::cout << " ( " << newDropped << " dropped this time)";
     441      }
     442      std::cout << std::endl;
     443    }
     444  }
     445  status_ = 0;
     446  return newDropped;
    452447}
    453448#else
    454449/* Factorize - filling in rowsDropped and returning number dropped */
    455 int
    456 ClpCholeskyUfl::factorize(const double * diagonal, int * rowsDropped)
     450int ClpCholeskyUfl::factorize(const double *diagonal, int *rowsDropped)
    457451{
    458452  return ClpCholeskyBase::factorize(diagonal, rowsDropped);
     
    462456#ifdef COIN_HAS_CHOLMOD
    463457/* Uses factorization to solve. */
    464 void
    465 ClpCholeskyUfl::solve (double * region)
    466 {
    467      cholmod_dense *x, *b;
    468      b = cholmod_allocate_dense (numberRows_, 1, numberRows_, CHOLMOD_REAL, c_) ;
    469      CoinMemcpyN(region, numberRows_, (double *) b->x);
    470      x = cholmod_solve (CHOLMOD_A, L_, b, c_) ;
    471      CoinMemcpyN((double *) x->x, numberRows_, region);
    472      cholmod_free_dense (&x, c_) ;
    473      cholmod_free_dense (&b, c_) ;
    474 }
    475 #else
    476 void
    477 ClpCholeskyUfl::solve (double * region)
     458void ClpCholeskyUfl::solve(double *region)
     459{
     460  cholmod_dense *x, *b;
     461  b = cholmod_allocate_dense(numberRows_, 1, numberRows_, CHOLMOD_REAL, c_);
     462  CoinMemcpyN(region, numberRows_, (double *)b->x);
     463  x = cholmod_solve(CHOLMOD_A, L_, b, c_);
     464  CoinMemcpyN((double *)x->x, numberRows_, region);
     465  cholmod_free_dense(&x, c_);
     466  cholmod_free_dense(&b, c_);
     467}
     468#else
     469void ClpCholeskyUfl::solve(double *region)
    478470{
    479471  ClpCholeskyBase::solve(region);
    480472}
    481473#endif
     474
     475/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     476*/
Note: See TracChangeset for help on using the changeset viewer.