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

    r1665 r2385  
    55// This code is licensed under the terms of the Eclipse Public License (EPL).
    66
    7 
    87#include "CoinPragma.hpp"
    98#include "CoinHelperFunctions.hpp"
     
    2120// Default Constructor
    2221//-------------------------------------------------------------------
    23 ClpCholeskyTaucs::ClpCholeskyTaucs ()
    24      : ClpCholeskyBase(),
    25        matrix_(NULL),
    26        factorization_(NULL),
    27        sparseFactorT_(NULL),
    28        choleskyStartT_(NULL),
    29        choleskyRowT_(NULL),
    30        sizeFactorT_(0),
    31       rowCopyT_(NULL)
    32 {
    33      type_ = 13;
     22ClpCholeskyTaucs::ClpCholeskyTaucs()
     23  : ClpCholeskyBase()
     24  , matrix_(NULL)
     25  , factorization_(NULL)
     26  , sparseFactorT_(NULL)
     27  , choleskyStartT_(NULL)
     28  , choleskyRowT_(NULL)
     29  , sizeFactorT_(0)
     30  , rowCopyT_(NULL)
     31{
     32  type_ = 13;
    3433}
    3534
     
    3736// Copy constructor
    3837//-------------------------------------------------------------------
    39 ClpCholeskyTaucs::ClpCholeskyTaucs (const ClpCholeskyTaucs & rhs)
    40      : ClpCholeskyBase(rhs)
    41 {
    42      type_ = rhs.type_;
    43      // For Taucs stuff is done by malloc
    44      matrix_ = rhs.matrix_;
    45      sizeFactorT_ = rhs.sizeFactorT_;
    46      if (matrix_) {
    47           choleskyStartT_ = (int *) malloc((numberRows_ + 1) * sizeof(int));
    48           CoinMemcpyN(rhs.choleskyStartT_, (numberRows_ + 1), choleskyStartT_);
    49           choleskyRowT_ = (int *) malloc(sizeFactorT_ * sizeof(int));
    50           CoinMemcpyN(rhs.choleskyRowT_, sizeFactorT_, choleskyRowT_);
    51           sparseFactorT_ = (double *) malloc(sizeFactorT_ * sizeof(double));
    52           CoinMemcpyN(rhs.sparseFactorT_, sizeFactorT_, sparseFactorT_);
    53           matrix_->colptr = choleskyStartT_;
    54           matrix_->rowind = choleskyRowT_;
    55           matrix_->values.d = sparseFactorT_;
    56      } else {
    57           sparseFactorT_ = NULL;
    58           choleskyStartT_ = NULL;
    59           choleskyRowT_ = NULL;
    60      }
    61      factorization_ = NULL,
    62      rowCopyT_ = rhs.rowCopyT_->clone();
    63 }
    64 
     38ClpCholeskyTaucs::ClpCholeskyTaucs(const ClpCholeskyTaucs &rhs)
     39  : ClpCholeskyBase(rhs)
     40{
     41  type_ = rhs.type_;
     42  // For Taucs stuff is done by malloc
     43  matrix_ = rhs.matrix_;
     44  sizeFactorT_ = rhs.sizeFactorT_;
     45  if (matrix_) {
     46    choleskyStartT_ = (int *)malloc((numberRows_ + 1) * sizeof(int));
     47    CoinMemcpyN(rhs.choleskyStartT_, (numberRows_ + 1), choleskyStartT_);
     48    choleskyRowT_ = (int *)malloc(sizeFactorT_ * sizeof(int));
     49    CoinMemcpyN(rhs.choleskyRowT_, sizeFactorT_, choleskyRowT_);
     50    sparseFactorT_ = (double *)malloc(sizeFactorT_ * sizeof(double));
     51    CoinMemcpyN(rhs.sparseFactorT_, sizeFactorT_, sparseFactorT_);
     52    matrix_->colptr = choleskyStartT_;
     53    matrix_->rowind = choleskyRowT_;
     54    matrix_->values.d = sparseFactorT_;
     55  } else {
     56    sparseFactorT_ = NULL;
     57    choleskyStartT_ = NULL;
     58    choleskyRowT_ = NULL;
     59  }
     60  factorization_ = NULL,
     61  rowCopyT_ = rhs.rowCopyT_->clone();
     62}
    6563
    6664//-------------------------------------------------------------------
    6765// Destructor
    6866//-------------------------------------------------------------------
    69 ClpCholeskyTaucs::~ClpCholeskyTaucs ()
    70 {
    71      taucs_ccs_free(matrix_);
    72      if (factorization_)
    73           taucs_supernodal_factor_free(factorization_);
    74      delete rowCopyT_;
     67ClpCholeskyTaucs::~ClpCholeskyTaucs()
     68{
     69  taucs_ccs_free(matrix_);
     70  if (factorization_)
     71    taucs_supernodal_factor_free(factorization_);
     72  delete rowCopyT_;
    7573}
    7674
     
    7977//-------------------------------------------------------------------
    8078ClpCholeskyTaucs &
    81 ClpCholeskyTaucs::operator=(const ClpCholeskyTaucs& rhs)
    82 {
    83      if (this != &rhs) {
    84           ClpCholeskyBase::operator=(rhs);
    85           taucs_ccs_free(matrix_);
    86           if (factorization_)
    87                taucs_supernodal_factor_free(factorization_);
    88           factorization_ = NULL;
    89           sizeFactorT_ = rhs.sizeFactorT_;
    90           matrix_ = rhs.matrix_;
    91           if (matrix_) {
    92                choleskyStartT_ = (int *) malloc((numberRows_ + 1) * sizeof(int));
    93                CoinMemcpyN(rhs.choleskyStartT_, (numberRows_ + 1), choleskyStartT_);
    94                choleskyRowT_ = (int *) malloc(sizeFactorT_ * sizeof(int));
    95                CoinMemcpyN(rhs.choleskyRowT_, sizeFactorT_, choleskyRowT_);
    96                sparseFactorT_ = (double *) malloc(sizeFactorT_ * sizeof(double));
    97                CoinMemcpyN(rhs.sparseFactorT_, sizeFactorT_, sparseFactorT_);
    98                matrix_->colptr = choleskyStartT_;
    99                matrix_->rowind = choleskyRowT_;
    100                matrix_->values.d = sparseFactorT_;
    101           } else {
    102                sparseFactorT_ = NULL;
    103                choleskyStartT_ = NULL;
    104                choleskyRowT_ = NULL;
     79ClpCholeskyTaucs::operator=(const ClpCholeskyTaucs &rhs)
     80{
     81  if (this != &rhs) {
     82    ClpCholeskyBase::operator=(rhs);
     83    taucs_ccs_free(matrix_);
     84    if (factorization_)
     85      taucs_supernodal_factor_free(factorization_);
     86    factorization_ = NULL;
     87    sizeFactorT_ = rhs.sizeFactorT_;
     88    matrix_ = rhs.matrix_;
     89    if (matrix_) {
     90      choleskyStartT_ = (int *)malloc((numberRows_ + 1) * sizeof(int));
     91      CoinMemcpyN(rhs.choleskyStartT_, (numberRows_ + 1), choleskyStartT_);
     92      choleskyRowT_ = (int *)malloc(sizeFactorT_ * sizeof(int));
     93      CoinMemcpyN(rhs.choleskyRowT_, sizeFactorT_, choleskyRowT_);
     94      sparseFactorT_ = (double *)malloc(sizeFactorT_ * sizeof(double));
     95      CoinMemcpyN(rhs.sparseFactorT_, sizeFactorT_, sparseFactorT_);
     96      matrix_->colptr = choleskyStartT_;
     97      matrix_->rowind = choleskyRowT_;
     98      matrix_->values.d = sparseFactorT_;
     99    } else {
     100      sparseFactorT_ = NULL;
     101      choleskyStartT_ = NULL;
     102      choleskyRowT_ = NULL;
     103    }
     104    delete rowCopyT_;
     105    rowCopyT_ = rhs.rowCopyT_->clone();
     106  }
     107  return *this;
     108}
     109//-------------------------------------------------------------------
     110// Clone
     111//-------------------------------------------------------------------
     112ClpCholeskyBase *ClpCholeskyTaucs::clone() const
     113{
     114  return new ClpCholeskyTaucs(*this);
     115}
     116/* Orders rows and saves pointer to matrix.and model */
     117int ClpCholeskyTaucs::order(ClpInterior *model)
     118{
     119  numberRows_ = model->numberRows();
     120  rowsDropped_ = new char[numberRows_];
     121  memset(rowsDropped_, 0, numberRows_);
     122  numberRowsDropped_ = 0;
     123  model_ = model;
     124  rowCopyT_ = model->clpMatrix()->reverseOrderedCopy();
     125  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
     126  const int *columnLength = model_->clpMatrix()->getVectorLengths();
     127  const int *row = model_->clpMatrix()->getIndices();
     128  const CoinBigIndex *rowStart = rowCopyT_->getVectorStarts();
     129  const int *rowLength = rowCopyT_->getVectorLengths();
     130  const int *column = rowCopyT_->getIndices();
     131  // We need two arrays for counts
     132  int *which = new int[numberRows_];
     133  int *used = new int[numberRows_];
     134  CoinZeroN(used, numberRows_);
     135  int iRow;
     136  sizeFactorT_ = 0;
     137  for (iRow = 0; iRow < numberRows_; iRow++) {
     138    int number = 1;
     139    // make sure diagonal exists
     140    which[0] = iRow;
     141    used[iRow] = 1;
     142    if (!rowsDropped_[iRow]) {
     143      CoinBigIndex startRow = rowStart[iRow];
     144      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
     145      for (CoinBigIndex k = startRow; k < endRow; k++) {
     146        int iColumn = column[k];
     147        CoinBigIndex start = columnStart[iColumn];
     148        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     149        for (CoinBigIndex j = start; j < end; j++) {
     150          int jRow = row[j];
     151          if (jRow >= iRow && !rowsDropped_[jRow]) {
     152            if (!used[jRow]) {
     153              used[jRow] = 1;
     154              which[number++] = jRow;
     155            }
    105156          }
    106           delete rowCopyT_;
    107           rowCopyT_ = rhs.rowCopyT_->clone();
    108      }
    109      return *this;
    110 }
    111 //-------------------------------------------------------------------
    112 // Clone
    113 //-------------------------------------------------------------------
    114 ClpCholeskyBase * ClpCholeskyTaucs::clone() const
    115 {
    116      return new ClpCholeskyTaucs(*this);
    117 }
    118 /* Orders rows and saves pointer to matrix.and model */
    119 int
    120 ClpCholeskyTaucs::order(ClpInterior * model)
    121 {
    122      numberRows_ = model->numberRows();
    123      rowsDropped_ = new char [numberRows_];
    124      memset(rowsDropped_, 0, numberRows_);
    125      numberRowsDropped_ = 0;
    126      model_ = model;
    127      rowCopyT_ = model->clpMatrix()->reverseOrderedCopy();
    128      const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    129      const int * columnLength = model_->clpMatrix()->getVectorLengths();
    130      const int * row = model_->clpMatrix()->getIndices();
    131      const CoinBigIndex * rowStart = rowCopyT_->getVectorStarts();
    132      const int * rowLength = rowCopyT_->getVectorLengths();
    133      const int * column = rowCopyT_->getIndices();
    134      // We need two arrays for counts
    135      int * which = new int [numberRows_];
    136      int * used = new int[numberRows_];
    137      CoinZeroN(used, numberRows_);
    138      int iRow;
    139      sizeFactorT_ = 0;
    140      for (iRow = 0; iRow < numberRows_; iRow++) {
    141           int number = 1;
    142           // make sure diagonal exists
    143           which[0] = iRow;
    144           used[iRow] = 1;
    145           if (!rowsDropped_[iRow]) {
    146                CoinBigIndex startRow = rowStart[iRow];
    147                CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
    148                for (CoinBigIndex k = startRow; k < endRow; k++) {
    149                     int iColumn = column[k];
    150                     CoinBigIndex start = columnStart[iColumn];
    151                     CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    152                     for (CoinBigIndex j = start; j < end; j++) {
    153                          int jRow = row[j];
    154                          if (jRow >= iRow && !rowsDropped_[jRow]) {
    155                               if (!used[jRow]) {
    156                                    used[jRow] = 1;
    157                                    which[number++] = jRow;
    158                               }
    159                          }
    160                     }
    161                }
    162                sizeFactorT_ += number;
    163                int j;
    164                for (j = 0; j < number; j++)
    165                     used[which[j]] = 0;
     157        }
     158      }
     159      sizeFactorT_ += number;
     160      int j;
     161      for (j = 0; j < number; j++)
     162        used[which[j]] = 0;
     163    }
     164  }
     165  delete[] which;
     166  // Now we have size - create arrays and fill in
     167  matrix_ = taucs_ccs_create(numberRows_, numberRows_, sizeFactorT_,
     168    TAUCS_DOUBLE | TAUCS_SYMMETRIC | TAUCS_LOWER);
     169  if (!matrix_)
     170    return 1;
     171  // Space for starts
     172  choleskyStartT_ = matrix_->colptr;
     173  choleskyRowT_ = matrix_->rowind;
     174  sparseFactorT_ = matrix_->values.d;
     175  sizeFactorT_ = 0;
     176  which = choleskyRowT_;
     177  for (iRow = 0; iRow < numberRows_; iRow++) {
     178    int number = 1;
     179    // make sure diagonal exists
     180    which[0] = iRow;
     181    used[iRow] = 1;
     182    choleskyStartT_[iRow] = sizeFactorT_;
     183    if (!rowsDropped_[iRow]) {
     184      CoinBigIndex startRow = rowStart[iRow];
     185      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
     186      for (CoinBigIndex k = startRow; k < endRow; k++) {
     187        int iColumn = column[k];
     188        CoinBigIndex start = columnStart[iColumn];
     189        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     190        for (CoinBigIndex j = start; j < end; j++) {
     191          int jRow = row[j];
     192          if (jRow >= iRow && !rowsDropped_[jRow]) {
     193            if (!used[jRow]) {
     194              used[jRow] = 1;
     195              which[number++] = jRow;
     196            }
    166197          }
    167      }
    168      delete [] which;
    169      // Now we have size - create arrays and fill in
    170      matrix_ = taucs_ccs_create(numberRows_, numberRows_, sizeFactorT_,
    171                                 TAUCS_DOUBLE | TAUCS_SYMMETRIC | TAUCS_LOWER);
    172      if (!matrix_)
    173           return 1;
    174      // Space for starts
    175      choleskyStartT_ = matrix_->colptr;
    176      choleskyRowT_ = matrix_->rowind;
    177      sparseFactorT_ = matrix_->values.d;
    178      sizeFactorT_ = 0;
    179      which = choleskyRowT_;
    180      for (iRow = 0; iRow < numberRows_; iRow++) {
    181           int number = 1;
    182           // make sure diagonal exists
    183           which[0] = iRow;
    184           used[iRow] = 1;
    185           choleskyStartT_[iRow] = sizeFactorT_;
    186           if (!rowsDropped_[iRow]) {
    187                CoinBigIndex startRow = rowStart[iRow];
    188                CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
    189                for (CoinBigIndex k = startRow; k < endRow; k++) {
    190                     int iColumn = column[k];
    191                     CoinBigIndex start = columnStart[iColumn];
    192                     CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    193                     for (CoinBigIndex j = start; j < end; j++) {
    194                          int jRow = row[j];
    195                          if (jRow >= iRow && !rowsDropped_[jRow]) {
    196                               if (!used[jRow]) {
    197                                    used[jRow] = 1;
    198                                    which[number++] = jRow;
    199                               }
    200                          }
    201                     }
    202                }
    203                sizeFactorT_ += number;
    204                int j;
    205                for (j = 0; j < number; j++)
    206                     used[which[j]] = 0;
    207                // Sort
    208                std::sort(which, which + number);
    209                // move which on
    210                which += number;
    211           }
    212      }
    213      choleskyStartT_[numberRows_] = sizeFactorT_;
    214      delete [] used;
    215      permuteInverse_ = new int [numberRows_];
    216      permute_ = new int[numberRows_];
    217      int * perm, *invp;
    218      // There seem to be bugs in ordering if model too small
    219      if (numberRows_ > 10)
    220           taucs_ccs_order(matrix_, &perm, &invp, (const char *) "genmmd");
    221      else
    222           taucs_ccs_order(matrix_, &perm, &invp, (const char *) "identity");
    223      CoinMemcpyN(perm, numberRows_, permuteInverse_);
    224      free(perm);
    225      CoinMemcpyN(invp, numberRows_, permute_);
    226      free(invp);
    227      // need to permute
    228      taucs_ccs_matrix * permuted = taucs_ccs_permute_symmetrically(matrix_, permuteInverse_, permute_);
    229      // symbolic
    230      factorization_ = taucs_ccs_factor_llt_symbolic(permuted);
    231      taucs_ccs_free(permuted);
    232      return factorization_ ? 0 : 1;
     198        }
     199      }
     200      sizeFactorT_ += number;
     201      int j;
     202      for (j = 0; j < number; j++)
     203        used[which[j]] = 0;
     204      // Sort
     205      std::sort(which, which + number);
     206      // move which on
     207      which += number;
     208    }
     209  }
     210  choleskyStartT_[numberRows_] = sizeFactorT_;
     211  delete[] used;
     212  permuteInverse_ = new int[numberRows_];
     213  permute_ = new int[numberRows_];
     214  int *perm, *invp;
     215  // There seem to be bugs in ordering if model too small
     216  if (numberRows_ > 10)
     217    taucs_ccs_order(matrix_, &perm, &invp, (const char *)"genmmd");
     218  else
     219    taucs_ccs_order(matrix_, &perm, &invp, (const char *)"identity");
     220  CoinMemcpyN(perm, numberRows_, permuteInverse_);
     221  free(perm);
     222  CoinMemcpyN(invp, numberRows_, permute_);
     223  free(invp);
     224  // need to permute
     225  taucs_ccs_matrix *permuted = taucs_ccs_permute_symmetrically(matrix_, permuteInverse_, permute_);
     226  // symbolic
     227  factorization_ = taucs_ccs_factor_llt_symbolic(permuted);
     228  taucs_ccs_free(permuted);
     229  return factorization_ ? 0 : 1;
    233230}
    234231/* Does Symbolic factorization given permutation.
     
    236233   user must provide factorize and solve.  Otherwise the default factorization is used
    237234   returns non-zero if not enough memory */
    238 int
    239 ClpCholeskyTaucs::symbolic()
    240 {
    241      return 0;
     235int ClpCholeskyTaucs::symbolic()
     236{
     237  return 0;
    242238}
    243239/* Factorize - filling in rowsDropped and returning number dropped */
    244 int
    245 ClpCholeskyTaucs::factorize(const double * diagonal, int * rowsDropped)
    246 {
    247      const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    248      const int * columnLength = model_->clpMatrix()->getVectorLengths();
    249      const int * row = model_->clpMatrix()->getIndices();
    250      const double * element = model_->clpMatrix()->getElements();
    251      const CoinBigIndex * rowStart = rowCopyT_->getVectorStarts();
    252      const int * rowLength = rowCopyT_->getVectorLengths();
    253      const int * column = rowCopyT_->getIndices();
    254      const double * elementByRow = rowCopyT_->getElements();
    255      int numberColumns = model_->clpMatrix()->getNumCols();
    256      int iRow;
    257      double * work = new double[numberRows_];
    258      CoinZeroN(work, numberRows_);
    259      const double * diagonalSlack = diagonal + numberColumns;
    260      int newDropped = 0;
    261      double largest;
    262      //perturbation
    263      double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
    264      perturbation = perturbation * perturbation;
    265      if (perturbation > 1.0) {
    266           //if (model_->model()->logLevel()&4)
    267           std::cout << "large perturbation " << perturbation << std::endl;
    268           perturbation = sqrt(perturbation);;
    269           perturbation = 1.0;
    270      }
    271      for (iRow = 0; iRow < numberRows_; iRow++) {
    272           double * put = sparseFactorT_ + choleskyStartT_[iRow];
    273           int * which = choleskyRowT_ + choleskyStartT_[iRow];
    274           int number = choleskyStartT_[iRow+1] - choleskyStartT_[iRow];
    275           if (!rowLength[iRow])
    276                rowsDropped_[iRow] = 1;
    277           if (!rowsDropped_[iRow]) {
    278                CoinBigIndex startRow = rowStart[iRow];
    279                CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
    280                work[iRow] = diagonalSlack[iRow];
    281                for (CoinBigIndex k = startRow; k < endRow; k++) {
    282                     int iColumn = column[k];
    283                     CoinBigIndex start = columnStart[iColumn];
    284                     CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    285                     double multiplier = diagonal[iColumn] * elementByRow[k];
    286                     for (CoinBigIndex j = start; j < end; j++) {
    287                          int jRow = row[j];
    288                          if (jRow >= iRow && !rowsDropped_[jRow]) {
    289                               double value = element[j] * multiplier;
    290                               work[jRow] += value;
    291                          }
    292                     }
    293                }
    294                int j;
    295                for (j = 0; j < number; j++) {
    296                     int jRow = which[j];
    297                     put[j] = work[jRow];
    298                     work[jRow] = 0.0;
    299                }
    300           } else {
    301                // dropped
    302                int j;
    303                for (j = 1; j < number; j++) {
    304                     put[j] = 0.0;
    305                }
    306                put[0] = 1.0;
     240int ClpCholeskyTaucs::factorize(const double *diagonal, int *rowsDropped)
     241{
     242  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
     243  const int *columnLength = model_->clpMatrix()->getVectorLengths();
     244  const int *row = model_->clpMatrix()->getIndices();
     245  const double *element = model_->clpMatrix()->getElements();
     246  const CoinBigIndex *rowStart = rowCopyT_->getVectorStarts();
     247  const int *rowLength = rowCopyT_->getVectorLengths();
     248  const int *column = rowCopyT_->getIndices();
     249  const double *elementByRow = rowCopyT_->getElements();
     250  int numberColumns = model_->clpMatrix()->getNumCols();
     251  int iRow;
     252  double *work = new double[numberRows_];
     253  CoinZeroN(work, numberRows_);
     254  const double *diagonalSlack = diagonal + numberColumns;
     255  int newDropped = 0;
     256  double largest;
     257  //perturbation
     258  double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
     259  perturbation = perturbation * perturbation;
     260  if (perturbation > 1.0) {
     261    //if (model_->model()->logLevel()&4)
     262    std::cout << "large perturbation " << perturbation << std::endl;
     263    perturbation = sqrt(perturbation);
     264    ;
     265    perturbation = 1.0;
     266  }
     267  for (iRow = 0; iRow < numberRows_; iRow++) {
     268    double *put = sparseFactorT_ + choleskyStartT_[iRow];
     269    int *which = choleskyRowT_ + choleskyStartT_[iRow];
     270    int number = choleskyStartT_[iRow + 1] - choleskyStartT_[iRow];
     271    if (!rowLength[iRow])
     272      rowsDropped_[iRow] = 1;
     273    if (!rowsDropped_[iRow]) {
     274      CoinBigIndex startRow = rowStart[iRow];
     275      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
     276      work[iRow] = diagonalSlack[iRow];
     277      for (CoinBigIndex k = startRow; k < endRow; k++) {
     278        int iColumn = column[k];
     279        CoinBigIndex start = columnStart[iColumn];
     280        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     281        double multiplier = diagonal[iColumn] * elementByRow[k];
     282        for (CoinBigIndex j = start; j < end; j++) {
     283          int jRow = row[j];
     284          if (jRow >= iRow && !rowsDropped_[jRow]) {
     285            double value = element[j] * multiplier;
     286            work[jRow] += value;
    307287          }
    308      }
    309      //check sizes
    310      double largest2 = maximumAbsElement(sparseFactorT_, sizeFactorT_);
    311      largest2 *= 1.0e-19;
    312      largest = CoinMin(largest2, 1.0e-11);
    313      int numberDroppedBefore = 0;
    314      for (iRow = 0; iRow < numberRows_; iRow++) {
    315           int dropped = rowsDropped_[iRow];
    316           // Move to int array
    317           rowsDropped[iRow] = dropped;
    318           if (!dropped) {
    319                CoinBigIndex start = choleskyStartT_[iRow];
    320                double diagonal = sparseFactorT_[start];
    321                if (diagonal > largest2) {
    322                     sparseFactorT_[start] = diagonal + perturbation;
    323                } else {
    324                     sparseFactorT_[start] = diagonal + perturbation;
    325                     rowsDropped[iRow] = 2;
    326                     numberDroppedBefore++;
    327                }
    328           }
    329      }
    330      taucs_supernodal_factor_free_numeric(factorization_);
    331      // need to permute
    332      taucs_ccs_matrix * permuted = taucs_ccs_permute_symmetrically(matrix_, permuteInverse_, permute_);
    333      int rCode = taucs_ccs_factor_llt_numeric(permuted, factorization_);
    334      taucs_ccs_free(permuted);
    335      if (rCode)
    336           printf("return code of %d from factor\n", rCode);
    337      delete [] work;
    338      choleskyCondition_ = 1.0;
    339      bool cleanCholesky;
    340      if (model_->numberIterations() < 200)
    341           cleanCholesky = true;
    342      else
    343           cleanCholesky = false;
    344      /*
     288        }
     289      }
     290      int j;
     291      for (j = 0; j < number; j++) {
     292        int jRow = which[j];
     293        put[j] = work[jRow];
     294        work[jRow] = 0.0;
     295      }
     296    } else {
     297      // dropped
     298      int j;
     299      for (j = 1; j < number; j++) {
     300        put[j] = 0.0;
     301      }
     302      put[0] = 1.0;
     303    }
     304  }
     305  //check sizes
     306  double largest2 = maximumAbsElement(sparseFactorT_, sizeFactorT_);
     307  largest2 *= 1.0e-19;
     308  largest = CoinMin(largest2, 1.0e-11);
     309  int numberDroppedBefore = 0;
     310  for (iRow = 0; iRow < numberRows_; iRow++) {
     311    int dropped = rowsDropped_[iRow];
     312    // Move to int array
     313    rowsDropped[iRow] = dropped;
     314    if (!dropped) {
     315      CoinBigIndex start = choleskyStartT_[iRow];
     316      double diagonal = sparseFactorT_[start];
     317      if (diagonal > largest2) {
     318        sparseFactorT_[start] = diagonal + perturbation;
     319      } else {
     320        sparseFactorT_[start] = diagonal + perturbation;
     321        rowsDropped[iRow] = 2;
     322        numberDroppedBefore++;
     323      }
     324    }
     325  }
     326  taucs_supernodal_factor_free_numeric(factorization_);
     327  // need to permute
     328  taucs_ccs_matrix *permuted = taucs_ccs_permute_symmetrically(matrix_, permuteInverse_, permute_);
     329  int rCode = taucs_ccs_factor_llt_numeric(permuted, factorization_);
     330  taucs_ccs_free(permuted);
     331  if (rCode)
     332    printf("return code of %d from factor\n", rCode);
     333  delete[] work;
     334  choleskyCondition_ = 1.0;
     335  bool cleanCholesky;
     336  if (model_->numberIterations() < 200)
     337    cleanCholesky = true;
     338  else
     339    cleanCholesky = false;
     340  /*
    345341       How do I find out where 1.0e100's are in cholesky?
    346342     */
    347      if (cleanCholesky) {
    348           //drop fresh makes some formADAT easier
    349           int oldDropped = numberRowsDropped_;
    350           if (newDropped || numberRowsDropped_) {
    351                std::cout << "Rank " << numberRows_ - newDropped << " ( " <<
    352                          newDropped << " dropped)";
    353                if (newDropped > oldDropped)
    354                     std::cout << " ( " << newDropped - oldDropped << " dropped this time)";
    355                std::cout << std::endl;
    356                newDropped = 0;
    357                for (int i = 0; i < numberRows_; i++) {
    358                     char dropped = rowsDropped[i];
    359                     rowsDropped_[i] = dropped;
    360                     if (dropped == 2) {
    361                          //dropped this time
    362                          rowsDropped[newDropped++] = i;
    363                          rowsDropped_[i] = 0;
    364                     }
    365                }
    366                numberRowsDropped_ = newDropped;
    367                newDropped = -(1 + newDropped);
    368           }
    369      } else {
    370           if (newDropped) {
    371                newDropped = 0;
    372                for (int i = 0; i < numberRows_; i++) {
    373                     char dropped = rowsDropped[i];
    374                     int oldDropped = rowsDropped_[i];;
    375                     rowsDropped_[i] = dropped;
    376                     if (dropped == 2) {
    377                          assert (!oldDropped);
    378                          //dropped this time
    379                          rowsDropped[newDropped++] = i;
    380                          rowsDropped_[i] = 1;
    381                     }
    382                }
    383           }
    384           numberRowsDropped_ += newDropped;
    385           if (numberRowsDropped_) {
    386                std::cout << "Rank " << numberRows_ - numberRowsDropped_ << " ( " <<
    387                          numberRowsDropped_ << " dropped)";
    388                if (newDropped) {
    389                     std::cout << " ( " << newDropped << " dropped this time)";
    390                }
    391                std::cout << std::endl;
    392           }
    393      }
    394      status_ = 0;
    395      return newDropped;
     343  if (cleanCholesky) {
     344    //drop fresh makes some formADAT easier
     345    int oldDropped = numberRowsDropped_;
     346    if (newDropped || numberRowsDropped_) {
     347      std::cout << "Rank " << numberRows_ - newDropped << " ( " << newDropped << " dropped)";
     348      if (newDropped > oldDropped)
     349        std::cout << " ( " << newDropped - oldDropped << " dropped this time)";
     350      std::cout << std::endl;
     351      newDropped = 0;
     352      for (int i = 0; i < numberRows_; i++) {
     353        char dropped = rowsDropped[i];
     354        rowsDropped_[i] = dropped;
     355        if (dropped == 2) {
     356          //dropped this time
     357          rowsDropped[newDropped++] = i;
     358          rowsDropped_[i] = 0;
     359        }
     360      }
     361      numberRowsDropped_ = newDropped;
     362      newDropped = -(1 + newDropped);
     363    }
     364  } else {
     365    if (newDropped) {
     366      newDropped = 0;
     367      for (int i = 0; i < numberRows_; i++) {
     368        char dropped = rowsDropped[i];
     369        int oldDropped = rowsDropped_[i];
     370        ;
     371        rowsDropped_[i] = dropped;
     372        if (dropped == 2) {
     373          assert(!oldDropped);
     374          //dropped this time
     375          rowsDropped[newDropped++] = i;
     376          rowsDropped_[i] = 1;
     377        }
     378      }
     379    }
     380    numberRowsDropped_ += newDropped;
     381    if (numberRowsDropped_) {
     382      std::cout << "Rank " << numberRows_ - numberRowsDropped_ << " ( " << numberRowsDropped_ << " dropped)";
     383      if (newDropped) {
     384        std::cout << " ( " << newDropped << " dropped this time)";
     385      }
     386      std::cout << std::endl;
     387    }
     388  }
     389  status_ = 0;
     390  return newDropped;
    396391}
    397392/* Uses factorization to solve. */
    398 void
    399 ClpCholeskyTaucs::solve (double * region)
    400 {
    401      double * in = new double[numberRows_];
    402      double * out = new double[numberRows_];
    403      taucs_vec_permute(numberRows_, TAUCS_DOUBLE, region, in, permuteInverse_);
    404      int rCode = taucs_supernodal_solve_llt(factorization_, out, in);
    405      if (rCode)
    406           printf("return code of %d from solve\n", rCode);
    407      taucs_vec_permute(numberRows_, TAUCS_DOUBLE, out, region, permute_);
    408      delete [] out;
    409      delete [] in;
     393void ClpCholeskyTaucs::solve(double *region)
     394{
     395  double *in = new double[numberRows_];
     396  double *out = new double[numberRows_];
     397  taucs_vec_permute(numberRows_, TAUCS_DOUBLE, region, in, permuteInverse_);
     398  int rCode = taucs_supernodal_solve_llt(factorization_, out, in);
     399  if (rCode)
     400    printf("return code of %d from solve\n", rCode);
     401  taucs_vec_permute(numberRows_, TAUCS_DOUBLE, out, region, permute_);
     402  delete[] out;
     403  delete[] in;
    410404}
    411405#endif
     406
     407/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     408*/
Note: See TracChangeset for help on using the changeset viewer.