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

    r2275 r2385  
    33// Corporation and others.  All Rights Reserved.
    44// This code is licensed under the terms of the Eclipse Public License (EPL).
    5 
    65
    76#include "CoinPragma.hpp"
     
    3635// Default Constructor
    3736//-------------------------------------------------------------------
    38 ClpCholeskyMumps::ClpCholeskyMumps (int denseThreshold,int logLevel)
    39      : ClpCholeskyBase(denseThreshold)
    40 {
    41      mumps_ = (DMUMPS_STRUC_C*)malloc(sizeof(DMUMPS_STRUC_C));
    42      type_ = 16;
    43      mumps_->n = 0;
    44      mumps_->nz = 0;
    45      mumps_->a = NULL;
    46      mumps_->jcn = NULL;
    47      mumps_->irn = NULL;
    48      mumps_->job = JOB_INIT;//initialize mumps
    49      mumps_->par = 1;//working host for sequential version
    50      mumps_->sym = 2;//general symmetric matrix
    51      mumps_->comm_fortran = USE_COMM_WORLD;
    52      int myid;
    53      int justName;
    54      MPI_Init(&justName, NULL);
     37ClpCholeskyMumps::ClpCholeskyMumps(int denseThreshold, int logLevel)
     38  : ClpCholeskyBase(denseThreshold)
     39{
     40  mumps_ = (DMUMPS_STRUC_C *)malloc(sizeof(DMUMPS_STRUC_C));
     41  type_ = 16;
     42  mumps_->n = 0;
     43  mumps_->nz = 0;
     44  mumps_->a = NULL;
     45  mumps_->jcn = NULL;
     46  mumps_->irn = NULL;
     47  mumps_->job = JOB_INIT; //initialize mumps
     48  mumps_->par = 1; //working host for sequential version
     49  mumps_->sym = 2; //general symmetric matrix
     50  mumps_->comm_fortran = USE_COMM_WORLD;
     51  int myid;
     52  int justName;
     53  MPI_Init(&justName, NULL);
    5554#ifndef NDEBUG
    56      int ierr = MPI_Comm_rank(MPI_COMM_WORLD, &myid);
    57      assert (!ierr);
     55  int ierr = MPI_Comm_rank(MPI_COMM_WORLD, &myid);
     56  assert(!ierr);
    5857#else
    59      MPI_Comm_rank(MPI_COMM_WORLD, &myid);
     58  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
    6059#endif
    61      dmumps_c(mumps_);
     60  dmumps_c(mumps_);
    6261#define ICNTL(I) icntl[(I)-1] /* macro s.t. indices match documentation */
    6362#define CNTL(I) cntl[(I)-1] /* macro s.t. indices match documentation */
    64      mumps_->ICNTL(5) = 1; // say compressed format
    65      mumps_->ICNTL(4) = 2; // log messages
    66      mumps_->ICNTL(24) = 1; // Deal with zeros on diagonal
    67      mumps_->CNTL(3) = 1.0e-20; // drop if diagonal less than this
    68      if (!logLevel) {
    69        // output off
    70        mumps_->ICNTL(1) = -1;
    71        mumps_->ICNTL(2) = -1;
    72        mumps_->ICNTL(3) = -1;
    73        mumps_->ICNTL(4) = 0;
    74      }
     63  mumps_->ICNTL(5) = 1; // say compressed format
     64  mumps_->ICNTL(4) = 2; // log messages
     65  mumps_->ICNTL(24) = 1; // Deal with zeros on diagonal
     66  mumps_->CNTL(3) = 1.0e-20; // drop if diagonal less than this
     67  if (!logLevel) {
     68    // output off
     69    mumps_->ICNTL(1) = -1;
     70    mumps_->ICNTL(2) = -1;
     71    mumps_->ICNTL(3) = -1;
     72    mumps_->ICNTL(4) = 0;
     73  }
    7574}
    7675
     
    7877// Copy constructor
    7978//-------------------------------------------------------------------
    80 ClpCholeskyMumps::ClpCholeskyMumps (const ClpCholeskyMumps & rhs)
    81      : ClpCholeskyBase(rhs)
    82 {
    83      abort();
    84 }
    85 
     79ClpCholeskyMumps::ClpCholeskyMumps(const ClpCholeskyMumps &rhs)
     80  : ClpCholeskyBase(rhs)
     81{
     82  abort();
     83}
    8684
    8785//-------------------------------------------------------------------
    8886// Destructor
    8987//-------------------------------------------------------------------
    90 ClpCholeskyMumps::~ClpCholeskyMumps ()
    91 {
    92      mumps_->job = JOB_END;
    93      dmumps_c(mumps_); /* Terminate instance */
    94      MPI_Finalize();
    95      free(mumps_);
     88ClpCholeskyMumps::~ClpCholeskyMumps()
     89{
     90  mumps_->job = JOB_END;
     91  dmumps_c(mumps_); /* Terminate instance */
     92  MPI_Finalize();
     93  free(mumps_);
    9694}
    9795
     
    10098//-------------------------------------------------------------------
    10199ClpCholeskyMumps &
    102 ClpCholeskyMumps::operator=(const ClpCholeskyMumps& rhs)
    103 {
    104      if (this != &rhs) {
    105           ClpCholeskyBase::operator=(rhs);
    106           abort();
    107      }
    108      return *this;
     100ClpCholeskyMumps::operator=(const ClpCholeskyMumps &rhs)
     101{
     102  if (this != &rhs) {
     103    ClpCholeskyBase::operator=(rhs);
     104    abort();
     105  }
     106  return *this;
    109107}
    110108//-------------------------------------------------------------------
    111109// Clone
    112110//-------------------------------------------------------------------
    113 ClpCholeskyBase * ClpCholeskyMumps::clone() const
    114 {
    115      return new ClpCholeskyMumps(*this);
     111ClpCholeskyBase *ClpCholeskyMumps::clone() const
     112{
     113  return new ClpCholeskyMumps(*this);
    116114}
    117115/* Orders rows and saves pointer to matrix.and model */
    118 int
    119 ClpCholeskyMumps::order(ClpInterior * model)
    120 {
    121      numberRows_ = model->numberRows();
    122      if (doKKT_) {
    123           numberRows_ += numberRows_ + model->numberColumns();
    124           printf("finish coding MUMPS KKT!\n");
    125           abort();
    126      }
    127      rowsDropped_ = new char [numberRows_];
    128      memset(rowsDropped_, 0, numberRows_);
    129      numberRowsDropped_ = 0;
    130      model_ = model;
    131      rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
    132      const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    133      const int * columnLength = model_->clpMatrix()->getVectorLengths();
    134      const int * row = model_->clpMatrix()->getIndices();
    135      const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
    136      const int * rowLength = rowCopy_->getVectorLengths();
    137      const int * column = rowCopy_->getIndices();
    138      // We need two arrays for counts
    139      int * which = new int [numberRows_];
    140      int * used = new int[numberRows_+1];
    141      CoinZeroN(used, numberRows_);
    142      int iRow;
    143      sizeFactor_ = 0;
    144      for (iRow = 0; iRow < numberRows_; iRow++) {
    145           int number = 1;
    146           // make sure diagonal exists
    147           which[0] = iRow;
    148           used[iRow] = 1;
    149           if (!rowsDropped_[iRow]) {
    150                CoinBigIndex startRow = rowStart[iRow];
    151                CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
    152                for (CoinBigIndex k = startRow; k < endRow; k++) {
    153                     int iColumn = column[k];
    154                     CoinBigIndex start = columnStart[iColumn];
    155                     CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    156                     for (CoinBigIndex j = start; j < end; j++) {
    157                          int jRow = row[j];
    158                          if (jRow >= iRow && !rowsDropped_[jRow]) {
    159                               if (!used[jRow]) {
    160                                    used[jRow] = 1;
    161                                    which[number++] = jRow;
    162                               }
    163                          }
    164                     }
    165                }
    166                sizeFactor_ += number;
    167                int j;
    168                for (j = 0; j < number; j++)
    169                     used[which[j]] = 0;
     116int ClpCholeskyMumps::order(ClpInterior *model)
     117{
     118  numberRows_ = model->numberRows();
     119  if (doKKT_) {
     120    numberRows_ += numberRows_ + model->numberColumns();
     121    printf("finish coding MUMPS KKT!\n");
     122    abort();
     123  }
     124  rowsDropped_ = new char[numberRows_];
     125  memset(rowsDropped_, 0, numberRows_);
     126  numberRowsDropped_ = 0;
     127  model_ = model;
     128  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
     129  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
     130  const int *columnLength = model_->clpMatrix()->getVectorLengths();
     131  const int *row = model_->clpMatrix()->getIndices();
     132  const CoinBigIndex *rowStart = rowCopy_->getVectorStarts();
     133  const int *rowLength = rowCopy_->getVectorLengths();
     134  const int *column = rowCopy_->getIndices();
     135  // We need two arrays for counts
     136  int *which = new int[numberRows_];
     137  int *used = new int[numberRows_ + 1];
     138  CoinZeroN(used, numberRows_);
     139  int iRow;
     140  sizeFactor_ = 0;
     141  for (iRow = 0; iRow < numberRows_; iRow++) {
     142    int number = 1;
     143    // make sure diagonal exists
     144    which[0] = iRow;
     145    used[iRow] = 1;
     146    if (!rowsDropped_[iRow]) {
     147      CoinBigIndex startRow = rowStart[iRow];
     148      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
     149      for (CoinBigIndex k = startRow; k < endRow; k++) {
     150        int iColumn = column[k];
     151        CoinBigIndex start = columnStart[iColumn];
     152        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     153        for (CoinBigIndex j = start; j < end; j++) {
     154          int jRow = row[j];
     155          if (jRow >= iRow && !rowsDropped_[jRow]) {
     156            if (!used[jRow]) {
     157              used[jRow] = 1;
     158              which[number++] = jRow;
     159            }
    170160          }
    171      }
    172      delete [] which;
    173      // NOT COMPRESSED FOR NOW ??? - Space for starts
    174      mumps_->ICNTL(5) = 0; // say NOT compressed format
    175      try {
    176           choleskyStart_ = new int[numberRows_+1+sizeFactor_];
    177      } catch (...) {
    178           // no memory
    179           return -1;
    180      }
    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                               }
    223                          }
    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;
     161        }
     162      }
     163      sizeFactor_ += number;
     164      int j;
     165      for (j = 0; j < number; j++)
     166        used[which[j]] = 0;
     167    }
     168  }
     169  delete[] which;
     170  // NOT COMPRESSED FOR NOW ??? - Space for starts
     171  mumps_->ICNTL(5) = 0; // say NOT compressed format
     172  try {
     173    choleskyStart_ = new int[numberRows_ + 1 + sizeFactor_];
     174  } catch (...) {
     175    // no memory
     176    return -1;
     177  }
     178  // Now we have size - create arrays and fill in
     179  try {
     180    choleskyRow_ = new int[sizeFactor_];
     181  } catch (...) {
     182    // no memory
     183    delete[] choleskyStart_;
     184    choleskyStart_ = NULL;
     185    return -1;
     186  }
     187  try {
     188    sparseFactor_ = new double[sizeFactor_];
     189  } catch (...) {
     190    // no memory
     191    delete[] choleskyRow_;
     192    choleskyRow_ = NULL;
     193    delete[] choleskyStart_;
     194    choleskyStart_ = NULL;
     195    return -1;
     196  }
     197
     198  sizeFactor_ = 0;
     199  which = choleskyRow_;
     200  for (iRow = 0; iRow < numberRows_; iRow++) {
     201    int number = 1;
     202    // make sure diagonal exists
     203    which[0] = iRow;
     204    used[iRow] = 1;
     205    choleskyStart_[iRow] = sizeFactor_;
     206    if (!rowsDropped_[iRow]) {
     207      CoinBigIndex startRow = rowStart[iRow];
     208      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
     209      for (CoinBigIndex k = startRow; k < endRow; k++) {
     210        int iColumn = column[k];
     211        CoinBigIndex start = columnStart[iColumn];
     212        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     213        for (CoinBigIndex j = start; j < end; j++) {
     214          int jRow = row[j];
     215          if (jRow >= iRow && !rowsDropped_[jRow]) {
     216            if (!used[jRow]) {
     217              used[jRow] = 1;
     218              which[number++] = jRow;
     219            }
    234220          }
    235      }
    236      choleskyStart_[numberRows_] = sizeFactor_;
    237      delete [] used;
    238      permuteInverse_ = new int [numberRows_];
    239      permute_ = new int[numberRows_];
    240      // To Fortran and fake
    241      for (iRow = 0; iRow < numberRows_ + 1; iRow++) {
    242           int k = choleskyStart_[iRow];
    243           int kEnd = choleskyStart_[iRow+1];
    244           k += numberRows_ + 1;
    245           kEnd += numberRows_ + 1;
    246           for (; k < kEnd; k++)
    247                choleskyStart_[k] = iRow + 1;
    248           choleskyStart_[iRow]++;
    249      }
    250      mumps_->nz = sizeFactor_;
    251      mumps_->irn = choleskyStart_ + numberRows_ + 1;
    252      mumps_->jcn = choleskyRow_;
    253      mumps_->a = NULL;
    254      for (CoinBigIndex i = 0; i < sizeFactor_; i++) {
    255           choleskyRow_[i]++;
     221        }
     222      }
     223      sizeFactor_ += number;
     224      int j;
     225      for (j = 0; j < number; j++)
     226        used[which[j]] = 0;
     227      // Sort
     228      std::sort(which, which + number);
     229      // move which on
     230      which += number;
     231    }
     232  }
     233  choleskyStart_[numberRows_] = sizeFactor_;
     234  delete[] used;
     235  permuteInverse_ = new int[numberRows_];
     236  permute_ = new int[numberRows_];
     237  // To Fortran and fake
     238  for (iRow = 0; iRow < numberRows_ + 1; iRow++) {
     239    int k = choleskyStart_[iRow];
     240    int kEnd = choleskyStart_[iRow + 1];
     241    k += numberRows_ + 1;
     242    kEnd += numberRows_ + 1;
     243    for (; k < kEnd; k++)
     244      choleskyStart_[k] = iRow + 1;
     245    choleskyStart_[iRow]++;
     246  }
     247  mumps_->nz = sizeFactor_;
     248  mumps_->irn = choleskyStart_ + numberRows_ + 1;
     249  mumps_->jcn = choleskyRow_;
     250  mumps_->a = NULL;
     251  for (CoinBigIndex i = 0; i < sizeFactor_; i++) {
     252    choleskyRow_[i]++;
    256253#ifndef NDEBUG
    257           assert (mumps_->irn[i] >= 1 && mumps_->irn[i] <= numberRows_);
    258           assert (mumps_->jcn[i] >= 1 && mumps_->jcn[i] <= numberRows_);
     254    assert(mumps_->irn[i] >= 1 && mumps_->irn[i] <= numberRows_);
     255    assert(mumps_->jcn[i] >= 1 && mumps_->jcn[i] <= numberRows_);
    259256#endif
    260      }
    261      // validate
    262      //mumps code here
    263      mumps_->n = numberRows_;
    264      mumps_->nelt = numberRows_;
    265      mumps_->eltptr = choleskyStart_;
    266      mumps_->eltvar = choleskyRow_;
    267      mumps_->a_elt = NULL;
    268      mumps_->rhs = NULL;
    269      mumps_->job = 1; // order
    270      dmumps_c(mumps_);
    271      mumps_->a = sparseFactor_;
    272      if (mumps_->infog[0]) {
    273           COIN_DETAIL_PRINT(printf("MUMPS ordering failed -error %d %d\n",
    274                                    mumps_->infog[0], mumps_->infog[1]));
    275           return 1;
    276      } else {
    277           double size = mumps_->infog[19];
    278           if (size < 0)
    279                size *= -1000000;
    280           COIN_DETAIL_PRINT(printf("%g nonzeros, flop count %g\n", size, mumps_->rinfog[0]));
    281      }
    282      for (iRow = 0; iRow < numberRows_; iRow++) {
    283           permuteInverse_[iRow] = iRow;
    284           permute_[iRow] = iRow;
    285      }
    286      return 0;
     257  }
     258  // validate
     259  //mumps code here
     260  mumps_->n = numberRows_;
     261  mumps_->nelt = numberRows_;
     262  mumps_->eltptr = choleskyStart_;
     263  mumps_->eltvar = choleskyRow_;
     264  mumps_->a_elt = NULL;
     265  mumps_->rhs = NULL;
     266  mumps_->job = 1; // order
     267  dmumps_c(mumps_);
     268  mumps_->a = sparseFactor_;
     269  if (mumps_->infog[0]) {
     270    COIN_DETAIL_PRINT(printf("MUMPS ordering failed -error %d %d\n",
     271      mumps_->infog[0], mumps_->infog[1]));
     272    return 1;
     273  } else {
     274    double size = mumps_->infog[19];
     275    if (size < 0)
     276      size *= -1000000;
     277    COIN_DETAIL_PRINT(printf("%g nonzeros, flop count %g\n", size, mumps_->rinfog[0]));
     278  }
     279  for (iRow = 0; iRow < numberRows_; iRow++) {
     280    permuteInverse_[iRow] = iRow;
     281    permute_[iRow] = iRow;
     282  }
     283  return 0;
    287284}
    288285/* Does Symbolic factorization given permutation.
     
    290287   user must provide factorize and solve.  Otherwise the default factorization is used
    291288   returns non-zero if not enough memory */
    292 int
    293 ClpCholeskyMumps::symbolic()
    294 {
    295      return 0;
     289int ClpCholeskyMumps::symbolic()
     290{
     291  return 0;
    296292}
    297293/* Factorize - filling in rowsDropped and returning number dropped */
    298 int
    299 ClpCholeskyMumps::factorize(const double * diagonal, int * rowsDropped)
    300 {
    301      const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    302      const int * columnLength = model_->clpMatrix()->getVectorLengths();
    303      const int * row = model_->clpMatrix()->getIndices();
    304      const double * element = model_->clpMatrix()->getElements();
    305      const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
    306      const int * rowLength = rowCopy_->getVectorLengths();
    307      const int * column = rowCopy_->getIndices();
    308      const double * elementByRow = rowCopy_->getElements();
    309      int numberColumns = model_->clpMatrix()->getNumCols();
    310      int iRow;
    311      double * work = new double[numberRows_];
    312      CoinZeroN(work, numberRows_);
    313      const double * diagonalSlack = diagonal + numberColumns;
    314      int newDropped = 0;
    315      //double smallest;
    316      //perturbation
    317      double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
    318      perturbation = 0.0;
    319      perturbation = perturbation * perturbation;
    320      if (perturbation > 1.0) {
     294int ClpCholeskyMumps::factorize(const double *diagonal, int *rowsDropped)
     295{
     296  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
     297  const int *columnLength = model_->clpMatrix()->getVectorLengths();
     298  const int *row = model_->clpMatrix()->getIndices();
     299  const double *element = model_->clpMatrix()->getElements();
     300  const CoinBigIndex *rowStart = rowCopy_->getVectorStarts();
     301  const int *rowLength = rowCopy_->getVectorLengths();
     302  const int *column = rowCopy_->getIndices();
     303  const double *elementByRow = rowCopy_->getElements();
     304  int numberColumns = model_->clpMatrix()->getNumCols();
     305  int iRow;
     306  double *work = new double[numberRows_];
     307  CoinZeroN(work, numberRows_);
     308  const double *diagonalSlack = diagonal + numberColumns;
     309  int newDropped = 0;
     310  //double smallest;
     311  //perturbation
     312  double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
     313  perturbation = 0.0;
     314  perturbation = perturbation * perturbation;
     315  if (perturbation > 1.0) {
    321316#ifdef COIN_DEVELOP
    322           //if (model_->model()->logLevel()&4)
    323           std::cout << "large perturbation " << perturbation << std::endl;
     317    //if (model_->model()->logLevel()&4)
     318    std::cout << "large perturbation " << perturbation << std::endl;
    324319#endif
    325           perturbation = sqrt(perturbation);;
    326           perturbation = 1.0;
    327      }
    328      double delta2 = model_->delta(); // add delta*delta to diagonal
    329      delta2 *= delta2;
    330      for (iRow = 0; iRow < numberRows_; iRow++) {
    331           double * put = sparseFactor_ + choleskyStart_[iRow] - 1; // Fortran
    332           int * which = choleskyRow_ + choleskyStart_[iRow] - 1; // Fortran
    333           int number = choleskyStart_[iRow+1] - choleskyStart_[iRow];
    334           if (!rowLength[iRow])
    335                rowsDropped_[iRow] = 1;
    336           if (!rowsDropped_[iRow]) {
    337                CoinBigIndex startRow = rowStart[iRow];
    338                CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
    339                work[iRow] = diagonalSlack[iRow] + delta2;
    340                for (CoinBigIndex k = startRow; k < endRow; k++) {
    341                     int iColumn = column[k];
    342                     if (!whichDense_ || !whichDense_[iColumn]) {
    343                          CoinBigIndex start = columnStart[iColumn];
    344                          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    345                          double multiplier = diagonal[iColumn] * elementByRow[k];
    346                          for (CoinBigIndex j = start; j < end; j++) {
    347                               int jRow = row[j];
    348                               if (jRow >= iRow && !rowsDropped_[jRow]) {
    349                                    double value = element[j] * multiplier;
    350                                    work[jRow] += value;
    351                               }
    352                          }
    353                     }
    354                }
    355                int j;
    356                for (j = 0; j < number; j++) {
    357                     int jRow = which[j] - 1; // to Fortran
    358                     put[j] = work[jRow];
    359                     work[jRow] = 0.0;
    360                }
    361           } else {
    362                // dropped
    363                int j;
    364                for (j = 1; j < number; j++) {
    365                     put[j] = 0.0;
    366                }
    367                put[0] = 1.0;
     320    perturbation = sqrt(perturbation);
     321    ;
     322    perturbation = 1.0;
     323  }
     324  double delta2 = model_->delta(); // add delta*delta to diagonal
     325  delta2 *= delta2;
     326  for (iRow = 0; iRow < numberRows_; iRow++) {
     327    double *put = sparseFactor_ + choleskyStart_[iRow] - 1; // Fortran
     328    int *which = choleskyRow_ + choleskyStart_[iRow] - 1; // Fortran
     329    int number = choleskyStart_[iRow + 1] - choleskyStart_[iRow];
     330    if (!rowLength[iRow])
     331      rowsDropped_[iRow] = 1;
     332    if (!rowsDropped_[iRow]) {
     333      CoinBigIndex startRow = rowStart[iRow];
     334      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
     335      work[iRow] = diagonalSlack[iRow] + delta2;
     336      for (CoinBigIndex k = startRow; k < endRow; k++) {
     337        int iColumn = column[k];
     338        if (!whichDense_ || !whichDense_[iColumn]) {
     339          CoinBigIndex start = columnStart[iColumn];
     340          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     341          double multiplier = diagonal[iColumn] * elementByRow[k];
     342          for (CoinBigIndex j = start; j < end; j++) {
     343            int jRow = row[j];
     344            if (jRow >= iRow && !rowsDropped_[jRow]) {
     345              double value = element[j] * multiplier;
     346              work[jRow] += value;
     347            }
    368348          }
    369      }
    370      //check sizes
    371      double largest2 = maximumAbsElement(sparseFactor_, sizeFactor_);
    372      largest2 *= 1.0e-20;
    373      int numberDroppedBefore = 0;
    374      for (iRow = 0; iRow < numberRows_; iRow++) {
    375           int dropped = rowsDropped_[iRow];
    376           // Move to int array
    377           rowsDropped[iRow] = dropped;
    378           if (!dropped) {
    379                int start = choleskyStart_[iRow] - 1; // to Fortran
    380                double diagonal = sparseFactor_[start];
    381                if (diagonal > largest2) {
    382                     sparseFactor_[start] = CoinMax(diagonal, 1.0e-10);
    383                } else {
    384                     sparseFactor_[start] = CoinMax(diagonal, 1.0e-10);
    385                     rowsDropped[iRow] = 2;
    386                     numberDroppedBefore++;
    387                }
    388           }
    389      }
    390      delete [] work;
    391      // code here
    392      mumps_->a_elt = sparseFactor_;
    393      mumps_->rhs = NULL;
    394      mumps_->job = 2; // factorize
    395      dmumps_c(mumps_);
    396      if (mumps_->infog[0]) {
    397           COIN_DETAIL_PRINT(printf("MUMPS factorization failed -error %d %d\n",
    398                                    mumps_->infog[0], mumps_->infog[1]));
    399      }
    400      choleskyCondition_ = 1.0;
    401      bool cleanCholesky;
    402      if (model_->numberIterations() < 2000)
    403           cleanCholesky = true;
    404      else
    405           cleanCholesky = false;
    406      if (cleanCholesky) {
    407           //drop fresh makes some formADAT easier
    408           //int oldDropped=numberRowsDropped_;
    409           if (newDropped || numberRowsDropped_) {
    410                //std::cout <<"Rank "<<numberRows_-newDropped<<" ( "<<
    411                //  newDropped<<" dropped)";
    412                //if (newDropped>oldDropped)
    413                //std::cout<<" ( "<<newDropped-oldDropped<<" dropped this time)";
    414                //std::cout<<std::endl;
    415                newDropped = 0;
    416                for (int i = 0; i < numberRows_; i++) {
    417                     int dropped = rowsDropped[i];
    418                     rowsDropped_[i] = (char)dropped;
    419                     if (dropped == 2) {
    420                          //dropped this time
    421                          rowsDropped[newDropped++] = i;
    422                          rowsDropped_[i] = 0;
    423                     }
    424                }
    425                numberRowsDropped_ = newDropped;
    426                newDropped = -(2 + newDropped);
    427           }
    428      } else {
    429           if (newDropped) {
    430                newDropped = 0;
    431                for (int i = 0; i < numberRows_; i++) {
    432                     int dropped = rowsDropped[i];
    433                     rowsDropped_[i] = (char)dropped;
    434                     if (dropped == 2) {
    435                          //dropped this time
    436                          rowsDropped[newDropped++] = i;
    437                          rowsDropped_[i] = 1;
    438                     }
    439                }
    440           }
    441           numberRowsDropped_ += newDropped;
    442           if (numberRowsDropped_ && 0) {
    443                std::cout << "Rank " << numberRows_ - numberRowsDropped_ << " ( " <<
    444                          numberRowsDropped_ << " dropped)";
    445                if (newDropped) {
    446                     std::cout << " ( " << newDropped << " dropped this time)";
    447                }
    448                std::cout << std::endl;
    449           }
    450      }
    451      status_ = 0;
    452      return newDropped;
     349        }
     350      }
     351      int j;
     352      for (j = 0; j < number; j++) {
     353        int jRow = which[j] - 1; // to Fortran
     354        put[j] = work[jRow];
     355        work[jRow] = 0.0;
     356      }
     357    } else {
     358      // dropped
     359      int j;
     360      for (j = 1; j < number; j++) {
     361        put[j] = 0.0;
     362      }
     363      put[0] = 1.0;
     364    }
     365  }
     366  //check sizes
     367  double largest2 = maximumAbsElement(sparseFactor_, sizeFactor_);
     368  largest2 *= 1.0e-20;
     369  int numberDroppedBefore = 0;
     370  for (iRow = 0; iRow < numberRows_; iRow++) {
     371    int dropped = rowsDropped_[iRow];
     372    // Move to int array
     373    rowsDropped[iRow] = dropped;
     374    if (!dropped) {
     375      int start = choleskyStart_[iRow] - 1; // to Fortran
     376      double diagonal = sparseFactor_[start];
     377      if (diagonal > largest2) {
     378        sparseFactor_[start] = CoinMax(diagonal, 1.0e-10);
     379      } else {
     380        sparseFactor_[start] = CoinMax(diagonal, 1.0e-10);
     381        rowsDropped[iRow] = 2;
     382        numberDroppedBefore++;
     383      }
     384    }
     385  }
     386  delete[] work;
     387  // code here
     388  mumps_->a_elt = sparseFactor_;
     389  mumps_->rhs = NULL;
     390  mumps_->job = 2; // factorize
     391  dmumps_c(mumps_);
     392  if (mumps_->infog[0]) {
     393    COIN_DETAIL_PRINT(printf("MUMPS factorization failed -error %d %d\n",
     394      mumps_->infog[0], mumps_->infog[1]));
     395  }
     396  choleskyCondition_ = 1.0;
     397  bool cleanCholesky;
     398  if (model_->numberIterations() < 2000)
     399    cleanCholesky = true;
     400  else
     401    cleanCholesky = false;
     402  if (cleanCholesky) {
     403    //drop fresh makes some formADAT easier
     404    //int oldDropped=numberRowsDropped_;
     405    if (newDropped || numberRowsDropped_) {
     406      //std::cout <<"Rank "<<numberRows_-newDropped<<" ( "<<
     407      //  newDropped<<" dropped)";
     408      //if (newDropped>oldDropped)
     409      //std::cout<<" ( "<<newDropped-oldDropped<<" dropped this time)";
     410      //std::cout<<std::endl;
     411      newDropped = 0;
     412      for (int i = 0; i < numberRows_; i++) {
     413        int dropped = rowsDropped[i];
     414        rowsDropped_[i] = (char)dropped;
     415        if (dropped == 2) {
     416          //dropped this time
     417          rowsDropped[newDropped++] = i;
     418          rowsDropped_[i] = 0;
     419        }
     420      }
     421      numberRowsDropped_ = newDropped;
     422      newDropped = -(2 + newDropped);
     423    }
     424  } else {
     425    if (newDropped) {
     426      newDropped = 0;
     427      for (int i = 0; i < numberRows_; i++) {
     428        int dropped = rowsDropped[i];
     429        rowsDropped_[i] = (char)dropped;
     430        if (dropped == 2) {
     431          //dropped this time
     432          rowsDropped[newDropped++] = i;
     433          rowsDropped_[i] = 1;
     434        }
     435      }
     436    }
     437    numberRowsDropped_ += newDropped;
     438    if (numberRowsDropped_ && 0) {
     439      std::cout << "Rank " << numberRows_ - numberRowsDropped_ << " ( " << numberRowsDropped_ << " dropped)";
     440      if (newDropped) {
     441        std::cout << " ( " << newDropped << " dropped this time)";
     442      }
     443      std::cout << std::endl;
     444    }
     445  }
     446  status_ = 0;
     447  return newDropped;
    453448}
    454449/* Uses factorization to solve. */
    455 void
    456 ClpCholeskyMumps::solve (double * region)
    457 {
    458      mumps_->rhs = region;
    459      mumps_->job = 3; // solve
    460      dmumps_c(mumps_);
    461 }
     450void ClpCholeskyMumps::solve(double *region)
     451{
     452  mumps_->rhs = region;
     453  mumps_->job = 3; // solve
     454  dmumps_c(mumps_);
     455}
     456
     457/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     458*/
Note: See TracChangeset for help on using the changeset viewer.