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

    r2236 r2385  
    2121// Default Constructor
    2222//-------------------------------------------------------------------
    23 ClpCholeskyPardiso::ClpCholeskyPardiso (int denseThreshold)
    24      : ClpCholeskyBase(denseThreshold)
    25 {
    26      type_ = 17;
     23ClpCholeskyPardiso::ClpCholeskyPardiso(int denseThreshold)
     24  : ClpCholeskyBase(denseThreshold)
     25{
     26  type_ = 17;
    2727}
    2828
     
    3030// Copy constructor
    3131//-------------------------------------------------------------------
    32 ClpCholeskyPardiso::ClpCholeskyPardiso (const ClpCholeskyPardiso & rhs)
    33      : ClpCholeskyBase(rhs)
    34 {
    35 }
    36 
     32ClpCholeskyPardiso::ClpCholeskyPardiso(const ClpCholeskyPardiso &rhs)
     33  : ClpCholeskyBase(rhs)
     34{
     35}
    3736
    3837//-------------------------------------------------------------------
    3938// Destructor
    4039//-------------------------------------------------------------------
    41 ClpCholeskyPardiso::~ClpCholeskyPardiso ()
     40ClpCholeskyPardiso::~ClpCholeskyPardiso()
    4241{
    4342  /* -------------------------------------------------------------------- */
    4443  /* .. Termination and release of memory. */
    4544  /* -------------------------------------------------------------------- */
    46   int phase = -1;           /* Release internal memory. */
    47   void * voidParameters=reinterpret_cast<void *>(doubleParameters_);
    48   MKL_INT mtype = -2;       /* Real symmetric matrix */
    49   MKL_INT nrhs = 1;     /* Number of right hand sides. */
    50   memset(integerParameters_,0,sizeof(integerParameters_));
    51   MKL_INT maxfct, mnum, error, msglvl=0;
     45  int phase = -1; /* Release internal memory. */
     46  void *voidParameters = reinterpret_cast< void * >(doubleParameters_);
     47  MKL_INT mtype = -2; /* Real symmetric matrix */
     48  MKL_INT nrhs = 1; /* Number of right hand sides. */
     49  memset(integerParameters_, 0, sizeof(integerParameters_));
     50  MKL_INT maxfct, mnum, error, msglvl = 0;
    5251  /* Auxiliary variables. */
    53   double ddum;          /* Double dummy */
    54   MKL_INT idum;         /* Integer dummy. */
    55   PARDISO (voidParameters, &maxfct, &mnum, &mtype, &phase,
    56            &numberRows_, sparseFactor_,
    57            choleskyStart_, choleskyRow_, &idum,
    58            &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error);
     52  double ddum; /* Double dummy */
     53  MKL_INT idum; /* Integer dummy. */
     54  PARDISO(voidParameters, &maxfct, &mnum, &mtype, &phase,
     55    &numberRows_, sparseFactor_,
     56    choleskyStart_, choleskyRow_, &idum,
     57    &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error);
    5958}
    6059
     
    6362//-------------------------------------------------------------------
    6463ClpCholeskyPardiso &
    65 ClpCholeskyPardiso::operator=(const ClpCholeskyPardiso& rhs)
    66 {
    67      if (this != &rhs) {
    68           ClpCholeskyBase::operator=(rhs);
    69      }
    70      return *this;
     64ClpCholeskyPardiso::operator=(const ClpCholeskyPardiso &rhs)
     65{
     66  if (this != &rhs) {
     67    ClpCholeskyBase::operator=(rhs);
     68  }
     69  return *this;
    7170}
    7271//-------------------------------------------------------------------
    7372// Clone
    7473//-------------------------------------------------------------------
    75 ClpCholeskyBase * ClpCholeskyPardiso::clone() const
    76 {
    77      return new ClpCholeskyPardiso(*this);
     74ClpCholeskyBase *ClpCholeskyPardiso::clone() const
     75{
     76  return new ClpCholeskyPardiso(*this);
    7877}
    7978/* Orders rows and saves pointer to matrix.and model */
    80 int
    81 ClpCholeskyPardiso::order(ClpInterior * model)
    82 {
    83      numberRows_ = model->numberRows();
    84      rowsDropped_ = new char [numberRows_];
    85      memset(rowsDropped_, 0, numberRows_);
    86      numberRowsDropped_ = 0;
    87      model_ = model;
    88      rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
    89      // Space for starts
    90      choleskyStart_ = new CoinBigIndex[numberRows_+1];
    91      const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    92      const int * columnLength = model_->clpMatrix()->getVectorLengths();
    93      const int * row = model_->clpMatrix()->getIndices();
    94      const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
    95      const int * rowLength = rowCopy_->getVectorLengths();
    96      const int * column = rowCopy_->getIndices();
    97      // We need two arrays for counts
    98      int * which = new int [numberRows_];
    99      int * used = new int[numberRows_+1];
    100      CoinZeroN(used, numberRows_);
    101      int iRow;
    102      sizeFactor_ = 0;
    103      int numberColumns = model->numberColumns();
    104      int numberDense = 0;
    105      denseThreshold_=0;
    106      if (denseThreshold_ > 0) {
    107           delete [] whichDense_;
    108           delete [] denseColumn_;
    109           delete dense_;
    110           whichDense_ = new char[numberColumns];
    111           int iColumn;
    112           used[numberRows_] = 0;
    113           for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    114                int length = columnLength[iColumn];
    115                used[length] += 1;
     79int ClpCholeskyPardiso::order(ClpInterior *model)
     80{
     81  numberRows_ = model->numberRows();
     82  rowsDropped_ = new char[numberRows_];
     83  memset(rowsDropped_, 0, numberRows_);
     84  numberRowsDropped_ = 0;
     85  model_ = model;
     86  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
     87  // Space for starts
     88  choleskyStart_ = new CoinBigIndex[numberRows_ + 1];
     89  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
     90  const int *columnLength = model_->clpMatrix()->getVectorLengths();
     91  const int *row = model_->clpMatrix()->getIndices();
     92  const CoinBigIndex *rowStart = rowCopy_->getVectorStarts();
     93  const int *rowLength = rowCopy_->getVectorLengths();
     94  const int *column = rowCopy_->getIndices();
     95  // We need two arrays for counts
     96  int *which = new int[numberRows_];
     97  int *used = new int[numberRows_ + 1];
     98  CoinZeroN(used, numberRows_);
     99  int iRow;
     100  sizeFactor_ = 0;
     101  int numberColumns = model->numberColumns();
     102  int numberDense = 0;
     103  denseThreshold_ = 0;
     104  if (denseThreshold_ > 0) {
     105    delete[] whichDense_;
     106    delete[] denseColumn_;
     107    delete dense_;
     108    whichDense_ = new char[numberColumns];
     109    int iColumn;
     110    used[numberRows_] = 0;
     111    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     112      int length = columnLength[iColumn];
     113      used[length] += 1;
     114    }
     115    int nLong = 0;
     116    int stop = CoinMax(denseThreshold_ / 2, 100);
     117    for (iRow = numberRows_; iRow >= stop; iRow--) {
     118      if (used[iRow])
     119        COIN_DETAIL_PRINT(printf("%d columns are of length %d\n", used[iRow], iRow));
     120      nLong += used[iRow];
     121      if (nLong > 50 || nLong > (numberColumns >> 2))
     122        break;
     123    }
     124    CoinZeroN(used, numberRows_);
     125    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     126      if (columnLength[iColumn] < denseThreshold_) {
     127        whichDense_[iColumn] = 0;
     128      } else {
     129        whichDense_[iColumn] = 1;
     130        numberDense++;
     131      }
     132    }
     133    if (!numberDense || numberDense > 100) {
     134      // free
     135      delete[] whichDense_;
     136      whichDense_ = NULL;
     137      denseColumn_ = NULL;
     138      dense_ = NULL;
     139    } else {
     140      // space for dense columns
     141      denseColumn_ = new double[numberDense * numberRows_];
     142      // dense cholesky
     143      dense_ = new ClpCholeskyDense();
     144      dense_->reserveSpace(NULL, numberDense);
     145      COIN_DETAIL_PRINT(printf("Taking %d columns as dense\n", numberDense));
     146    }
     147  }
     148  for (iRow = 0; iRow < numberRows_; iRow++) {
     149    int number = 1;
     150    // make sure diagonal exists
     151    which[0] = iRow;
     152    used[iRow] = 1;
     153    if (!rowsDropped_[iRow]) {
     154      CoinBigIndex startRow = rowStart[iRow];
     155      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
     156      for (CoinBigIndex k = startRow; k < endRow; k++) {
     157        int iColumn = column[k];
     158        if (!whichDense_ || !whichDense_[iColumn]) {
     159          CoinBigIndex start = columnStart[iColumn];
     160          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     161          for (CoinBigIndex j = start; j < end; j++) {
     162            int jRow = row[j];
     163            if (jRow >= iRow && !rowsDropped_[jRow]) {
     164              if (!used[jRow]) {
     165                used[jRow] = 1;
     166                which[number++] = jRow;
     167              }
     168            }
    116169          }
    117           int nLong = 0;
    118           int stop = CoinMax(denseThreshold_ / 2, 100);
    119           for (iRow = numberRows_; iRow >= stop; iRow--) {
    120                if (used[iRow])
    121                  COIN_DETAIL_PRINT(printf("%d columns are of length %d\n", used[iRow], iRow));
    122                nLong += used[iRow];
    123                if (nLong > 50 || nLong > (numberColumns >> 2))
    124                     break;
     170        }
     171      }
     172      sizeFactor_ += number;
     173      int j;
     174      for (j = 0; j < number; j++)
     175        used[which[j]] = 0;
     176    }
     177  }
     178  delete[] which;
     179  // Now we have size - create arrays and fill in
     180  try {
     181    choleskyRow_ = new int[sizeFactor_];
     182  } catch (...) {
     183    // no memory
     184    delete[] choleskyStart_;
     185    choleskyStart_ = NULL;
     186    return -1;
     187  }
     188  try {
     189    sparseFactor_ = new double[sizeFactor_];
     190  } catch (...) {
     191    // no memory
     192    delete[] choleskyRow_;
     193    choleskyRow_ = NULL;
     194    delete[] choleskyStart_;
     195    choleskyStart_ = NULL;
     196    return -1;
     197  }
     198
     199  sizeFactor_ = 0;
     200  which = choleskyRow_;
     201  for (iRow = 0; iRow < numberRows_; iRow++) {
     202    int number = 1;
     203    // make sure diagonal exists
     204    which[0] = iRow;
     205    used[iRow] = 1;
     206    choleskyStart_[iRow] = sizeFactor_;
     207    if (!rowsDropped_[iRow]) {
     208      CoinBigIndex startRow = rowStart[iRow];
     209      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
     210      for (CoinBigIndex k = startRow; k < endRow; k++) {
     211        int iColumn = column[k];
     212        if (!whichDense_ || !whichDense_[iColumn]) {
     213          CoinBigIndex start = columnStart[iColumn];
     214          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     215          for (CoinBigIndex j = start; j < end; j++) {
     216            int jRow = row[j];
     217            if (jRow >= iRow && !rowsDropped_[jRow]) {
     218              if (!used[jRow]) {
     219                used[jRow] = 1;
     220                which[number++] = jRow;
     221              }
     222            }
    125223          }
    126           CoinZeroN(used, numberRows_);
    127           for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    128                if (columnLength[iColumn] < denseThreshold_) {
    129                     whichDense_[iColumn] = 0;
    130                } else {
    131                     whichDense_[iColumn] = 1;
    132                     numberDense++;
    133                }
    134           }
    135           if (!numberDense || numberDense > 100) {
    136                // free
    137                delete [] whichDense_;
    138                whichDense_ = NULL;
    139                denseColumn_ = NULL;
    140                dense_ = NULL;
    141           } else {
    142                // space for dense columns
    143                denseColumn_ = new double [numberDense*numberRows_];
    144                // dense cholesky
    145                dense_ = new ClpCholeskyDense();
    146                dense_->reserveSpace(NULL, numberDense);
    147                COIN_DETAIL_PRINT(printf("Taking %d columns as dense\n", numberDense));
    148           }
    149      }
    150      for (iRow = 0; iRow < numberRows_; iRow++) {
    151           int number = 1;
    152           // make sure diagonal exists
    153           which[0] = iRow;
    154           used[iRow] = 1;
    155           if (!rowsDropped_[iRow]) {
    156                CoinBigIndex startRow = rowStart[iRow];
    157                CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
    158                for (CoinBigIndex k = startRow; k < endRow; k++) {
    159                     int iColumn = column[k];
    160                     if (!whichDense_ || !whichDense_[iColumn]) {
    161                          CoinBigIndex start = columnStart[iColumn];
    162                          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    163                          for (CoinBigIndex j = start; j < end; j++) {
    164                               int jRow = row[j];
    165                               if (jRow >= iRow && !rowsDropped_[jRow]) {
    166                                    if (!used[jRow]) {
    167                                         used[jRow] = 1;
    168                                         which[number++] = jRow;
    169                                    }
    170                               }
    171                          }
    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                     if (!whichDense_ || !whichDense_[iColumn]) {
    215                          CoinBigIndex start = columnStart[iColumn];
    216                          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    217                          for (CoinBigIndex j = start; j < end; j++) {
    218                               int jRow = row[j];
    219                               if (jRow >= iRow && !rowsDropped_[jRow]) {
    220                                    if (!used[jRow]) {
    221                                         used[jRow] = 1;
    222                                         which[number++] = jRow;
    223                                    }
    224                               }
    225                          }
    226                     }
    227                }
    228                sizeFactor_ += number;
    229                int j;
    230                for (j = 0; j < number; j++)
    231                     used[which[j]] = 0;
    232                // Sort
    233                std::sort(which, which + number);
    234                // move which on
    235                which += number;
    236           }
    237      }
    238      choleskyStart_[numberRows_] = sizeFactor_;
    239      delete [] used;
    240      permuteInverse_ = new int [numberRows_];
    241      permute_ = new int[numberRows_];
    242      // Assume void * same size as double
    243      void * voidParameters=reinterpret_cast<void *>(doubleParameters_);
    244      MKL_INT mtype = -2;       /* Real symmetric matrix */
    245      MKL_INT nrhs = 1;     /* Number of right hand sides. */
    246      memset(integerParameters_,0,sizeof(integerParameters_));
    247      MKL_INT maxfct, mnum, phase, error, msglvl;
    248      /* Auxiliary variables. */
    249      double ddum;          /* Double dummy */
    250      MKL_INT idum;         /* Integer dummy. */
    251      /* -------------------------------------------------------------------- */
    252      /* .. Setup Pardiso control parameters. */
    253      /* -------------------------------------------------------------------- */
    254     integerParameters_[0] = 1;         /* No solver default */
    255     integerParameters_[1] = 2;         /* Fill-in reordering from METIS */
    256     integerParameters_[3] = 0;         /* No iterative-direct algorithm */
    257     integerParameters_[4] = 0;         /* No user fill-in reducing permutation */
    258     integerParameters_[5] = 0;         /* Write solution into x */
    259     integerParameters_[6] = 0;         /* Number of refinements done */
    260     integerParameters_[7] = 2;         /* Max numbers of iterative refinement steps */
    261     integerParameters_[8] = 0;         /* Not in use */
    262     integerParameters_[9] = 13;        /* Perturb the pivot elements with 1E-13 */
    263     integerParameters_[10] = 1;        /* Use nonsymmetric permutation and scaling MPS */
    264     integerParameters_[11] = 0;        /* Not in use */
    265     integerParameters_[12] = 1;        /* Maximum weighted matching algorithm is switched-off (default for symmetric). Try integerParameters_[12] = 1 in case of inappropriate accuracy */
    266     integerParameters_[13] = 0;        /* Output: Number of perturbed pivots */
    267     integerParameters_[14] = 0;        /* Not in use */
    268     integerParameters_[15] = 0;        /* Not in use */
    269     integerParameters_[16] = 0;        /* Not in use */
    270     integerParameters_[17] = -1;       /* Output: Number of nonzeros in the factor LU */
    271     integerParameters_[18] = -1;       /* Output: Mflops for LU factorization */
    272     integerParameters_[19] = 0;        /* Output: Numbers of CG Iterations */
    273     integerParameters_[19] = 1;        // for 2x2 pivoting
    274     integerParameters_[34] = 1;        /* PARDISO use C-style indexing for ia and ja arrays */
    275     integerParameters_[55] = 1;        /* Pivoting control */
    276     maxfct = 1;           /* Maximum number of numerical factorizations. */
    277     mnum = 1;         /* Which factorization to use. */
    278     msglvl = 0;           /* Print statistical information in file */
    279     error = 0;            /* Initialize error flag */
    280 /* -------------------------------------------------------------------- */
    281 /* .. Initialize the internal solver memory pointer. This is only */
    282 /* necessary for the FIRST call of the PARDISO solver. */
    283 /* -------------------------------------------------------------------- */
    284      memset(voidParameters,0,sizeof(doubleParameters_));
    285 /* -------------------------------------------------------------------- */
    286 /* .. Reordering and Symbolic Factorization. This step also allocates */
    287 /* all memory that is necessary for the factorization. */
    288 /* -------------------------------------------------------------------- */
    289      // temp - fill elements
    290      for (int i=0;i<sizeFactor_;i++)
    291        sparseFactor_[i]=1.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  // Assume void * same size as double
     241  void *voidParameters = reinterpret_cast< void * >(doubleParameters_);
     242  MKL_INT mtype = -2; /* Real symmetric matrix */
     243  MKL_INT nrhs = 1; /* Number of right hand sides. */
     244  memset(integerParameters_, 0, sizeof(integerParameters_));
     245  MKL_INT maxfct, mnum, phase, error, msglvl;
     246  /* Auxiliary variables. */
     247  double ddum; /* Double dummy */
     248  MKL_INT idum; /* Integer dummy. */
     249  /* -------------------------------------------------------------------- */
     250  /* .. Setup Pardiso control parameters. */
     251  /* -------------------------------------------------------------------- */
     252  integerParameters_[0] = 1; /* No solver default */
     253  integerParameters_[1] = 2; /* Fill-in reordering from METIS */
     254  integerParameters_[3] = 0; /* No iterative-direct algorithm */
     255  integerParameters_[4] = 0; /* No user fill-in reducing permutation */
     256  integerParameters_[5] = 0; /* Write solution into x */
     257  integerParameters_[6] = 0; /* Number of refinements done */
     258  integerParameters_[7] = 2; /* Max numbers of iterative refinement steps */
     259  integerParameters_[8] = 0; /* Not in use */
     260  integerParameters_[9] = 13; /* Perturb the pivot elements with 1E-13 */
     261  integerParameters_[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
     262  integerParameters_[11] = 0; /* Not in use */
     263  integerParameters_[12] = 1; /* Maximum weighted matching algorithm is switched-off (default for symmetric). Try integerParameters_[12] = 1 in case of inappropriate accuracy */
     264  integerParameters_[13] = 0; /* Output: Number of perturbed pivots */
     265  integerParameters_[14] = 0; /* Not in use */
     266  integerParameters_[15] = 0; /* Not in use */
     267  integerParameters_[16] = 0; /* Not in use */
     268  integerParameters_[17] = -1; /* Output: Number of nonzeros in the factor LU */
     269  integerParameters_[18] = -1; /* Output: Mflops for LU factorization */
     270  integerParameters_[19] = 0; /* Output: Numbers of CG Iterations */
     271  integerParameters_[19] = 1; // for 2x2 pivoting
     272  integerParameters_[34] = 1; /* PARDISO use C-style indexing for ia and ja arrays */
     273  integerParameters_[55] = 1; /* Pivoting control */
     274  maxfct = 1; /* Maximum number of numerical factorizations. */
     275  mnum = 1; /* Which factorization to use. */
     276  msglvl = 0; /* Print statistical information in file */
     277  error = 0; /* Initialize error flag */
     278  /* -------------------------------------------------------------------- */
     279  /* .. Initialize the internal solver memory pointer. This is only */
     280  /* necessary for the FIRST call of the PARDISO solver. */
     281  /* -------------------------------------------------------------------- */
     282  memset(voidParameters, 0, sizeof(doubleParameters_));
     283  /* -------------------------------------------------------------------- */
     284  /* .. Reordering and Symbolic Factorization. This step also allocates */
     285  /* all memory that is necessary for the factorization. */
     286  /* -------------------------------------------------------------------- */
     287  // temp - fill elements
     288  for (int i = 0; i < sizeFactor_; i++)
     289    sparseFactor_[i] = 1.0;
    292290#define ALWAYS_ORDER
    293291#ifdef ALWAYS_ORDER
    294      error=0;
     292  error = 0;
    295293#else
    296     phase = 11;
    297     PARDISO (voidParameters, &maxfct, &mnum, &mtype, &phase,
    298              &numberRows_, sparseFactor_,
    299              choleskyStart_, choleskyRow_, &idum,
    300              &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error);
     294  phase = 11;
     295  PARDISO(voidParameters, &maxfct, &mnum, &mtype, &phase,
     296    &numberRows_, sparseFactor_,
     297    choleskyStart_, choleskyRow_, &idum,
     298    &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error);
    301299#endif
    302     if ( error != 0 )
    303     {
    304         printf ("ERROR during symbolic factorization: %d\n", error);
    305         return 1;
    306     }
    307     printf ("Number of nonzeros in factors = %d\n", integerParameters_[17]);
    308     printf ("Number of factorization MFLOPS = %d\n", integerParameters_[18]);
    309     lastNumberDropped_=0;
    310      // Modify gamma and delta if 0.0
    311      if (!model->gamma() && !model->delta()) {
    312           model->setGamma(5.0e-5);
    313           model->setDelta(5.0e-5);
    314      }
    315      return 0;
     300  if (error != 0) {
     301    printf("ERROR during symbolic factorization: %d\n", error);
     302    return 1;
     303  }
     304  printf("Number of nonzeros in factors = %d\n", integerParameters_[17]);
     305  printf("Number of factorization MFLOPS = %d\n", integerParameters_[18]);
     306  lastNumberDropped_ = 0;
     307  // Modify gamma and delta if 0.0
     308  if (!model->gamma() && !model->delta()) {
     309    model->setGamma(5.0e-5);
     310    model->setDelta(5.0e-5);
     311  }
     312  return 0;
    316313}
    317314/* Does Symbolic factorization given permutation.
     
    319316   user must provide factorize and solve.  Otherwise the default factorization is used
    320317   returns non-zero if not enough memory */
    321 int
    322 ClpCholeskyPardiso::symbolic()
    323 {
    324      return 0;
     318int ClpCholeskyPardiso::symbolic()
     319{
     320  return 0;
    325321}
    326322/* Factorize - filling in rowsDropped and returning number dropped */
    327 int
    328 ClpCholeskyPardiso::factorize(const double * diagonal, int * rowsDropped)
    329 {
    330      const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    331      const int * columnLength = model_->clpMatrix()->getVectorLengths();
    332      const int * row = model_->clpMatrix()->getIndices();
    333      const double * element = model_->clpMatrix()->getElements();
    334      const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
    335      const int * rowLength = rowCopy_->getVectorLengths();
    336      const int * column = rowCopy_->getIndices();
    337      const double * elementByRow = rowCopy_->getElements();
    338      int numberColumns = model_->clpMatrix()->getNumCols();
    339      int iRow;
    340      double * work = new double[numberRows_];
    341      CoinZeroN(work, numberRows_);
    342      const double * diagonalSlack = diagonal + numberColumns;
    343      double largest;
    344      //int numberDense = 0;
    345      //if (dense_)
    346      //   numberDense = dense_->numberRows();
    347      //perturbation
    348      double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
    349      perturbation = perturbation * perturbation;
    350      if (perturbation > 1.0) {
     323int ClpCholeskyPardiso::factorize(const double *diagonal, int *rowsDropped)
     324{
     325  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
     326  const int *columnLength = model_->clpMatrix()->getVectorLengths();
     327  const int *row = model_->clpMatrix()->getIndices();
     328  const double *element = model_->clpMatrix()->getElements();
     329  const CoinBigIndex *rowStart = rowCopy_->getVectorStarts();
     330  const int *rowLength = rowCopy_->getVectorLengths();
     331  const int *column = rowCopy_->getIndices();
     332  const double *elementByRow = rowCopy_->getElements();
     333  int numberColumns = model_->clpMatrix()->getNumCols();
     334  int iRow;
     335  double *work = new double[numberRows_];
     336  CoinZeroN(work, numberRows_);
     337  const double *diagonalSlack = diagonal + numberColumns;
     338  double largest;
     339  //int numberDense = 0;
     340  //if (dense_)
     341  //   numberDense = dense_->numberRows();
     342  //perturbation
     343  double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
     344  perturbation = perturbation * perturbation;
     345  if (perturbation > 1.0) {
    351346#ifdef COIN_DEVELOP
    352           //if (model_->model()->logLevel()&4)
    353           std::cout << "large perturbation " << perturbation << std::endl;
     347    //if (model_->model()->logLevel()&4)
     348    std::cout << "large perturbation " << perturbation << std::endl;
    354349#endif
    355           perturbation = sqrt(perturbation);;
    356           perturbation = 1.0;
    357      }
    358      if (whichDense_) {
    359           double * denseDiagonal = dense_->diagonal();
    360           double * dense = denseColumn_;
    361           int iDense = 0;
    362           for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    363                if (whichDense_[iColumn]) {
    364                     denseDiagonal[iDense++] = 1.0 / diagonal[iColumn];
    365                     CoinZeroN(dense, numberRows_);
    366                     CoinBigIndex start = columnStart[iColumn];
    367                     CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    368                     for (CoinBigIndex j = start; j < end; j++) {
    369                          int jRow = row[j];
    370                          dense[jRow] = element[j];
    371                     }
    372                     dense += numberRows_;
    373                }
     350    perturbation = sqrt(perturbation);
     351    ;
     352    perturbation = 1.0;
     353  }
     354  if (whichDense_) {
     355    double *denseDiagonal = dense_->diagonal();
     356    double *dense = denseColumn_;
     357    int iDense = 0;
     358    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     359      if (whichDense_[iColumn]) {
     360        denseDiagonal[iDense++] = 1.0 / diagonal[iColumn];
     361        CoinZeroN(dense, numberRows_);
     362        CoinBigIndex start = columnStart[iColumn];
     363        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     364        for (CoinBigIndex j = start; j < end; j++) {
     365          int jRow = row[j];
     366          dense[jRow] = element[j];
     367        }
     368        dense += numberRows_;
     369      }
     370    }
     371  }
     372  double delta2 = model_->delta(); // add delta*delta to diagonal
     373  delta2 *= delta2;
     374  int *whichNon = new int[numberRows_];
     375  int nNon = 0;
     376  for (iRow = 0; iRow < numberRows_; iRow++) {
     377    //for (int i=0;i<numberRows_;i++)
     378    //assert(!work[i]);
     379    for (int i = 0; i < nNon; i++)
     380      work[whichNon[i]] = 0.0;
     381    nNon = 1;
     382    whichNon[0] = iRow;
     383    double *put = sparseFactor_ + choleskyStart_[iRow];
     384    int *which = choleskyRow_ + choleskyStart_[iRow];
     385    int number = choleskyStart_[iRow + 1] - choleskyStart_[iRow];
     386    if (!rowLength[iRow])
     387      rowsDropped_[iRow] = 1;
     388    if (!rowsDropped_[iRow]) {
     389      CoinBigIndex startRow = rowStart[iRow];
     390      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
     391      work[iRow] = diagonalSlack[iRow] + delta2;
     392      for (CoinBigIndex k = startRow; k < endRow; k++) {
     393        int iColumn = column[k];
     394        if (!whichDense_ || !whichDense_[iColumn]) {
     395          CoinBigIndex start = columnStart[iColumn];
     396          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     397          double multiplier = diagonal[iColumn] * elementByRow[k];
     398          for (CoinBigIndex j = start; j < end; j++) {
     399            int jRow = row[j];
     400            if (jRow >= iRow && !rowsDropped_[jRow]) {
     401              double value = element[j] * multiplier;
     402              if (!work[jRow])
     403                whichNon[nNon++] = jRow;
     404              work[jRow] += value;
     405            }
    374406          }
    375      }
    376      double delta2 = model_->delta(); // add delta*delta to diagonal
    377      delta2 *= delta2;
    378      int * whichNon = new int [numberRows_];
    379      int nNon=0;
    380      for (iRow = 0; iRow < numberRows_; iRow++) {
    381        //for (int i=0;i<numberRows_;i++)
    382        //assert(!work[i]);
    383        for (int i=0;i<nNon;i++)
    384          work[whichNon[i]]=0.0;
    385        nNon=1;
    386        whichNon[0]=iRow;
    387           double * put = sparseFactor_ + choleskyStart_[iRow];
    388           int * which = choleskyRow_ + choleskyStart_[iRow];
    389           int number = choleskyStart_[iRow+1] - choleskyStart_[iRow];
    390           if (!rowLength[iRow])
    391                rowsDropped_[iRow] = 1;
    392           if (!rowsDropped_[iRow]) {
    393                CoinBigIndex startRow = rowStart[iRow];
    394                CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
    395                work[iRow] = diagonalSlack[iRow] + delta2;
    396                for (CoinBigIndex k = startRow; k < endRow; k++) {
    397                     int iColumn = column[k];
    398                     if (!whichDense_ || !whichDense_[iColumn]) {
    399                          CoinBigIndex start = columnStart[iColumn];
    400                          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    401                          double multiplier = diagonal[iColumn] * elementByRow[k];
    402                          for (CoinBigIndex j = start; j < end; j++) {
    403                               int jRow = row[j];
    404                               if (jRow >= iRow && !rowsDropped_[jRow]) {
    405                                    double value = element[j] * multiplier;
    406                                    if (!work[jRow])
    407                                      whichNon[nNon++]=jRow;
    408                                    work[jRow] += value;
    409                               }
    410                          }
    411                     }
    412                }
    413                int j;
    414                for (j = 0; j < number; j++) {
    415                     int jRow = which[j];
    416                     put[j] = work[jRow];
    417                     work[jRow] = 0.0;
    418                }
    419           } else {
    420                // dropped
    421                int j;
    422                for (j = 1; j < number; j++) {
    423                     put[j] = 0.0;
    424                }
    425                put[0] = 1.0e12;
     407        }
     408      }
     409      int j;
     410      for (j = 0; j < number; j++) {
     411        int jRow = which[j];
     412        put[j] = work[jRow];
     413        work[jRow] = 0.0;
     414      }
     415    } else {
     416      // dropped
     417      int j;
     418      for (j = 1; j < number; j++) {
     419        put[j] = 0.0;
     420      }
     421      put[0] = 1.0e12;
     422    }
     423  }
     424  delete[] whichNon;
     425  //check sizes
     426  double largest2 = maximumAbsElement(sparseFactor_, sizeFactor_);
     427  largest2 *= 1.0e-20;
     428  largest = CoinMin(largest2, 1.0e-11);
     429  int numberDroppedBefore = 0;
     430  for (iRow = 0; iRow < numberRows_; iRow++) {
     431    int dropped = rowsDropped_[iRow];
     432    // Move to int array
     433    rowsDropped[iRow] = dropped;
     434    if (!dropped) {
     435      CoinBigIndex start = choleskyStart_[iRow];
     436      double diagonal = sparseFactor_[start];
     437      if (diagonal > largest2) {
     438        sparseFactor_[start] = diagonal + perturbation;
     439      } else {
     440        sparseFactor_[start] = diagonal + perturbation;
     441        rowsDropped[iRow] = 2;
     442        numberDroppedBefore++;
     443      }
     444    }
     445  }
     446  //integerParameters_[0] = 0;
     447  int maxfct = 1; /* Maximum number of numerical factorizations. */
     448  int mnum = 1; /* Which factorization to use. */
     449  int msglvl = 0; /* Print statistical information in file */
     450  int error = 0; /* Initialize error flag */
     451  int phase = 22;
     452#ifdef ALWAYS_ORDER
     453  {
     454    int phase = -1; /* Release internal memory. */
     455    void *voidParameters = reinterpret_cast< void * >(doubleParameters_);
     456    MKL_INT mtype = -2; /* Real symmetric matrix */
     457    MKL_INT nrhs = 1; /* Number of right hand sides. */
     458    memset(integerParameters_, 0, sizeof(integerParameters_));
     459    MKL_INT maxfct, mnum, error, msglvl = 0;
     460    /* Auxiliary variables. */
     461    double ddum; /* Double dummy */
     462    MKL_INT idum; /* Integer dummy. */
     463    PARDISO(voidParameters, &maxfct, &mnum, &mtype, &phase,
     464      &numberRows_, sparseFactor_,
     465      choleskyStart_, choleskyRow_, &idum,
     466      &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error);
     467    memset(integerParameters_, 0, sizeof(integerParameters_));
     468    /* -------------------------------------------------------------------- */
     469    /* .. Setup Pardiso control parameters. */
     470    /* -------------------------------------------------------------------- */
     471    integerParameters_[0] = 1; /* No solver default */
     472    integerParameters_[1] = 2; /* Fill-in reordering from METIS */
     473    integerParameters_[3] = 0; /* No iterative-direct algorithm */
     474    integerParameters_[4] = 0; /* No user fill-in reducing permutation */
     475    integerParameters_[5] = 0; /* Write solution into x */
     476    integerParameters_[6] = 0; /* Number of refinements done */
     477    integerParameters_[7] = 2; /* Max numbers of iterative refinement steps */
     478    integerParameters_[8] = 0; /* Not in use */
     479    integerParameters_[9] = 13; /* Perturb the pivot elements with 1E-13 */
     480    integerParameters_[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
     481    integerParameters_[11] = 0; /* Not in use */
     482    integerParameters_[12] = 1; /* Maximum weighted matching algorithm is switched-off (default for symmetric). Try integerParameters_[12] = 1 in case of inappropriate accuracy */
     483    integerParameters_[13] = 0; /* Output: Number of perturbed pivots */
     484    integerParameters_[14] = 0; /* Not in use */
     485    integerParameters_[15] = 0; /* Not in use */
     486    integerParameters_[16] = 0; /* Not in use */
     487    integerParameters_[17] = -1; /* Output: Number of nonzeros in the factor LU */
     488    integerParameters_[18] = -1; /* Output: Mflops for LU factorization */
     489    integerParameters_[19] = 0; /* Output: Numbers of CG Iterations */
     490    integerParameters_[19] = 1; // for 2x2 pivoting
     491    integerParameters_[34] = 1; /* PARDISO use C-style indexing for ia and ja arrays */
     492    integerParameters_[55] = 1; /* Pivoting control */
     493    maxfct = 1; /* Maximum number of numerical factorizations. */
     494    mnum = 1; /* Which factorization to use. */
     495    msglvl = 0; /* Print statistical information in file */
     496    error = 0; /* Initialize error flag */
     497    /* -------------------------------------------------------------------- */
     498    /* .. Initialize the internal solver memory pointer. This is only */
     499    /* necessary for the FIRST call of the PARDISO solver. */
     500    /* -------------------------------------------------------------------- */
     501    memset(voidParameters, 0, sizeof(doubleParameters_));
     502    /* -------------------------------------------------------------------- */
     503    /* .. Reordering and Symbolic Factorization. This step also allocates */
     504    /* all memory that is necessary for the factorization. */
     505    /* -------------------------------------------------------------------- */
     506    // pack down
     507    CoinBigIndex numberElements = 0;
     508    for (int iRow = 0; iRow < numberRows_; iRow++) {
     509      CoinBigIndex start = choleskyStart_[iRow];
     510      CoinBigIndex end = choleskyStart_[iRow + 1];
     511      choleskyStart_[iRow] = numberElements;
     512      for (CoinBigIndex j = start; j < end; j++) {
     513        if (sparseFactor_[j]) {
     514          choleskyRow_[numberElements] = choleskyRow_[j];
     515          sparseFactor_[numberElements++] = sparseFactor_[j];
     516        }
     517      }
     518    }
     519    choleskyStart_[numberRows_] = numberElements;
     520    phase = 11;
     521    PARDISO(voidParameters, &maxfct, &mnum, &mtype, &phase,
     522      &numberRows_, sparseFactor_,
     523      choleskyStart_, choleskyRow_, &idum,
     524      &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error);
     525    if (error != 0) {
     526      printf("ERROR during symbolic factorization: %d\n", error);
     527      return 1;
     528    }
     529    printf("Number of nonzeros in factors = %d\n", integerParameters_[17]);
     530    printf("Number of factorization MFLOPS = %d\n", integerParameters_[18]);
     531  }
     532#endif
     533  if (numberRowsDropped_ > lastNumberDropped_ && false) {
     534    phase = 12; //reorder
     535    CoinBigIndex numberElements = 0;
     536    for (int iRow = 0; iRow < numberRows_; iRow++) {
     537      CoinBigIndex start = choleskyStart_[iRow];
     538      CoinBigIndex end = choleskyStart_[iRow + 1];
     539      choleskyStart_[iRow] = numberElements;
     540      if (rowsDropped_[iRow] != 2) {
     541        for (CoinBigIndex j = start; j < end; j++) {
     542          choleskyRow_[numberElements] = choleskyRow_[j];
     543          sparseFactor_[numberElements++] = sparseFactor_[j];
     544        }
     545      } else {
     546        rowsDropped_[iRow] = 1;
     547        for (CoinBigIndex j = start; j < end; j++) {
     548          if (sparseFactor_[j]) {
     549            choleskyRow_[numberElements] = choleskyRow_[j];
     550            sparseFactor_[numberElements++] = sparseFactor_[j];
    426551          }
    427      }
    428      delete [] whichNon;
    429      //check sizes
    430      double largest2 = maximumAbsElement(sparseFactor_, sizeFactor_);
    431      largest2 *= 1.0e-20;
    432      largest = CoinMin(largest2, 1.0e-11);
    433      int numberDroppedBefore = 0;
    434      for (iRow = 0; iRow < numberRows_; iRow++) {
    435           int dropped = rowsDropped_[iRow];
    436           // Move to int array
    437           rowsDropped[iRow] = dropped;
    438           if (!dropped) {
    439                CoinBigIndex start = choleskyStart_[iRow];
    440                double diagonal = sparseFactor_[start];
    441                if (diagonal > largest2) {
    442                     sparseFactor_[start] = diagonal + perturbation;
    443                } else {
    444                     sparseFactor_[start] = diagonal + perturbation;
    445                     rowsDropped[iRow] = 2;
    446                     numberDroppedBefore++;
    447                }
    448           }
    449      }
    450      //integerParameters_[0] = 0;
    451      int maxfct = 1;           /* Maximum number of numerical factorizations. */
    452      int mnum = 1;         /* Which factorization to use. */
    453      int msglvl = 0;           /* Print statistical information in file */
    454      int error = 0;            /* Initialize error flag */
    455      int phase = 22;
    456 #ifdef ALWAYS_ORDER
    457      {
    458        int phase = -1;           /* Release internal memory. */
    459        void * voidParameters=reinterpret_cast<void *>(doubleParameters_);
    460        MKL_INT mtype = -2;       /* Real symmetric matrix */
    461        MKL_INT nrhs = 1;     /* Number of right hand sides. */
    462        memset(integerParameters_,0,sizeof(integerParameters_));
    463        MKL_INT maxfct, mnum, error, msglvl=0;
    464        /* Auxiliary variables. */
    465        double ddum;          /* Double dummy */
    466        MKL_INT idum;         /* Integer dummy. */
    467        PARDISO (voidParameters, &maxfct, &mnum, &mtype, &phase,
    468                 &numberRows_, sparseFactor_,
    469                 choleskyStart_, choleskyRow_, &idum,
    470                 &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error);
    471        memset(integerParameters_,0,sizeof(integerParameters_));
    472        /* -------------------------------------------------------------------- */
    473        /* .. Setup Pardiso control parameters. */
    474        /* -------------------------------------------------------------------- */
    475        integerParameters_[0] = 1;         /* No solver default */
    476        integerParameters_[1] = 2;         /* Fill-in reordering from METIS */
    477        integerParameters_[3] = 0;         /* No iterative-direct algorithm */
    478        integerParameters_[4] = 0;         /* No user fill-in reducing permutation */
    479        integerParameters_[5] = 0;         /* Write solution into x */
    480        integerParameters_[6] = 0;         /* Number of refinements done */
    481        integerParameters_[7] = 2;         /* Max numbers of iterative refinement steps */
    482        integerParameters_[8] = 0;         /* Not in use */
    483        integerParameters_[9] = 13;        /* Perturb the pivot elements with 1E-13 */
    484        integerParameters_[10] = 1;        /* Use nonsymmetric permutation and scaling MPS */
    485        integerParameters_[11] = 0;        /* Not in use */
    486        integerParameters_[12] = 1;        /* Maximum weighted matching algorithm is switched-off (default for symmetric). Try integerParameters_[12] = 1 in case of inappropriate accuracy */
    487        integerParameters_[13] = 0;        /* Output: Number of perturbed pivots */
    488        integerParameters_[14] = 0;        /* Not in use */
    489        integerParameters_[15] = 0;        /* Not in use */
    490        integerParameters_[16] = 0;        /* Not in use */
    491        integerParameters_[17] = -1;       /* Output: Number of nonzeros in the factor LU */
    492        integerParameters_[18] = -1;       /* Output: Mflops for LU factorization */
    493        integerParameters_[19] = 0;        /* Output: Numbers of CG Iterations */
    494        integerParameters_[19] = 1;        // for 2x2 pivoting
    495        integerParameters_[34] = 1;        /* PARDISO use C-style indexing for ia and ja arrays */
    496        integerParameters_[55] = 1;        /* Pivoting control */
    497        maxfct = 1;           /* Maximum number of numerical factorizations. */
    498        mnum = 1;         /* Which factorization to use. */
    499        msglvl = 0;           /* Print statistical information in file */
    500        error = 0;            /* Initialize error flag */
    501        /* -------------------------------------------------------------------- */
    502        /* .. Initialize the internal solver memory pointer. This is only */
    503        /* necessary for the FIRST call of the PARDISO solver. */
    504        /* -------------------------------------------------------------------- */
    505        memset(voidParameters,0,sizeof(doubleParameters_));
    506        /* -------------------------------------------------------------------- */
    507        /* .. Reordering and Symbolic Factorization. This step also allocates */
    508        /* all memory that is necessary for the factorization. */
    509        /* -------------------------------------------------------------------- */
    510        // pack down
    511        CoinBigIndex numberElements=0;
    512        for (int iRow = 0; iRow < numberRows_; iRow++) {
    513          CoinBigIndex start = choleskyStart_[iRow];
    514          CoinBigIndex end = choleskyStart_[iRow+1];
    515          choleskyStart_[iRow]=numberElements;
    516          for (CoinBigIndex j=start;j<end;j++) {
    517            if (sparseFactor_[j]) {
    518              choleskyRow_[numberElements]=choleskyRow_[j];
    519              sparseFactor_[numberElements++]=sparseFactor_[j];
    520            }
    521          }
    522        }
    523        choleskyStart_[numberRows_]=numberElements;
    524        phase = 11;
    525        PARDISO (voidParameters, &maxfct, &mnum, &mtype, &phase,
    526                 &numberRows_, sparseFactor_,
    527                 choleskyStart_, choleskyRow_, &idum,
    528              &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error);
    529        if ( error != 0 )
    530          {
    531            printf ("ERROR during symbolic factorization: %d\n", error);
    532            return 1;
    533          }
    534        printf ("Number of nonzeros in factors = %d\n", integerParameters_[17]);
    535        printf ("Number of factorization MFLOPS = %d\n", integerParameters_[18]);
    536      }
    537 #endif
    538      if (numberRowsDropped_>lastNumberDropped_&&false) {
    539        phase=12; //reorder
    540        CoinBigIndex numberElements=0;
    541        for (int iRow = 0; iRow < numberRows_; iRow++) {
    542          CoinBigIndex start = choleskyStart_[iRow];
    543          CoinBigIndex end = choleskyStart_[iRow+1];
    544          choleskyStart_[iRow]=numberElements;
    545          if (rowsDropped_[iRow]!=2) {
    546            for (CoinBigIndex j=start;j<end;j++) {
    547              choleskyRow_[numberElements]=choleskyRow_[j];
    548              sparseFactor_[numberElements++]=sparseFactor_[j];
    549            }
    550          } else {
    551            rowsDropped_[iRow]=1;
    552            for (CoinBigIndex j=start;j<end;j++) {
    553              if (sparseFactor_[j]) {
    554                choleskyRow_[numberElements]=choleskyRow_[j];
    555                sparseFactor_[numberElements++]=sparseFactor_[j];
    556              }
    557            }
    558          }
    559        }
    560        choleskyStart_[numberRows_]=numberElements;
    561        lastNumberDropped_=numberRowsDropped_;
    562      }
    563      MKL_INT mtype = -2;       /* Real symmetric matrix */
    564      void * voidParameters=reinterpret_cast<void *>(doubleParameters_);
    565      MKL_INT nrhs = 1;     /* Number of right hand sides. */
    566      /* Auxiliary variables. */
    567      double ddum;          /* Double dummy */
    568      MKL_INT idum;         /* Integer dummy. */
    569      PARDISO (voidParameters, &maxfct, &mnum, &mtype, &phase,
    570              &numberRows_, sparseFactor_,
    571              choleskyStart_, choleskyRow_, &idum,
    572              &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error);
    573      printf("%d positive, %d negative\n",integerParameters_[21],
    574             integerParameters_[22]);
    575      // use some other array?
    576      double * originalDiagonal=new double[numberRows_];
    577      pardiso_getdiag(voidParameters,work,originalDiagonal,&mnum,&error);
    578      delete [] originalDiagonal;
    579      largest=1.0e-11;
    580      double smallest=COIN_DBL_MAX;
    581      int newDropped = 0;
    582      for (int i=0;i<numberRows_;i++) {
    583        double value=work[i];
    584        if (value!=1.0e100) {
    585          //if (value>1.0e13)
    586          //printf("large diagonal %g at %d\n",value,i);
    587          largest=CoinMax(largest,value);
    588          smallest=CoinMin(smallest,value);
    589        } else if (!rowsDropped_[i]) {
    590          rowsDropped_[i]=2;
    591          rowsDropped[i]=2;
    592          printf("%d dropped\n",i);
    593          newDropped++;
    594          numberRowsDropped_++;
    595        }
    596      }
    597      if(smallest<=0.0)
    598        printf("small negative diagonal %g\n",smallest);
    599      delete [] work;
    600      if (model_->messageHandler()->logLevel() > 1)
    601           std::cout << "Cholesky - largest " << largest << " smallest " << smallest << std::endl;
    602      choleskyCondition_ = largest / smallest;
    603      status_ = 0;
    604      return newDropped;
     552        }
     553      }
     554    }
     555    choleskyStart_[numberRows_] = numberElements;
     556    lastNumberDropped_ = numberRowsDropped_;
     557  }
     558  MKL_INT mtype = -2; /* Real symmetric matrix */
     559  void *voidParameters = reinterpret_cast< void * >(doubleParameters_);
     560  MKL_INT nrhs = 1; /* Number of right hand sides. */
     561  /* Auxiliary variables. */
     562  double ddum; /* Double dummy */
     563  MKL_INT idum; /* Integer dummy. */
     564  PARDISO(voidParameters, &maxfct, &mnum, &mtype, &phase,
     565    &numberRows_, sparseFactor_,
     566    choleskyStart_, choleskyRow_, &idum,
     567    &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error);
     568  printf("%d positive, %d negative\n", integerParameters_[21],
     569    integerParameters_[22]);
     570  // use some other array?
     571  double *originalDiagonal = new double[numberRows_];
     572  pardiso_getdiag(voidParameters, work, originalDiagonal, &mnum, &error);
     573  delete[] originalDiagonal;
     574  largest = 1.0e-11;
     575  double smallest = COIN_DBL_MAX;
     576  int newDropped = 0;
     577  for (int i = 0; i < numberRows_; i++) {
     578    double value = work[i];
     579    if (value != 1.0e100) {
     580      //if (value>1.0e13)
     581      //printf("large diagonal %g at %d\n",value,i);
     582      largest = CoinMax(largest, value);
     583      smallest = CoinMin(smallest, value);
     584    } else if (!rowsDropped_[i]) {
     585      rowsDropped_[i] = 2;
     586      rowsDropped[i] = 2;
     587      printf("%d dropped\n", i);
     588      newDropped++;
     589      numberRowsDropped_++;
     590    }
     591  }
     592  if (smallest <= 0.0)
     593    printf("small negative diagonal %g\n", smallest);
     594  delete[] work;
     595  if (model_->messageHandler()->logLevel() > 1)
     596    std::cout << "Cholesky - largest " << largest << " smallest " << smallest << std::endl;
     597  choleskyCondition_ = largest / smallest;
     598  status_ = 0;
     599  return newDropped;
    605600}
    606601/* Uses factorization to solve. */
    607 void
    608 ClpCholeskyPardiso::solve (double * region)
     602void ClpCholeskyPardiso::solve(double *region)
    609603{
    610604  /* -------------------------------------------------------------------- */
     
    612606  /* -------------------------------------------------------------------- */
    613607  int phase = 33;
    614   integerParameters_[7] = 2;         /* Max numbers of iterative refinement steps. */
    615   double * x = new double[numberRows_];
    616   memcpy(x,region,numberRows_*sizeof(double));
    617   MKL_INT mtype = -2;       /* Real symmetric matrix */
    618   void * voidParameters=reinterpret_cast<void *>(doubleParameters_);
    619   MKL_INT nrhs = 1;     /* Number of right hand sides. */
    620   MKL_INT maxfct=1, mnum=1, error, msglvl=0;
     608  integerParameters_[7] = 2; /* Max numbers of iterative refinement steps. */
     609  double *x = new double[numberRows_];
     610  memcpy(x, region, numberRows_ * sizeof(double));
     611  MKL_INT mtype = -2; /* Real symmetric matrix */
     612  void *voidParameters = reinterpret_cast< void * >(doubleParameters_);
     613  MKL_INT nrhs = 1; /* Number of right hand sides. */
     614  MKL_INT maxfct = 1, mnum = 1, error, msglvl = 0;
    621615  /* Auxiliary variables. */
    622   MKL_INT idum;         /* Integer dummy. */
    623   PARDISO (voidParameters, &maxfct, &mnum, &mtype, &phase,
    624            &numberRows_, sparseFactor_,
    625            choleskyStart_, choleskyRow_, &idum,
    626            &nrhs, integerParameters_, &msglvl, x,region, &error);
    627   if ( error != 0 )
    628     {
    629       printf ("ERROR during solution: %d\n", error);
    630       exit (3);
    631     }
    632   delete [] x;
     616  MKL_INT idum; /* Integer dummy. */
     617  PARDISO(voidParameters, &maxfct, &mnum, &mtype, &phase,
     618    &numberRows_, sparseFactor_,
     619    choleskyStart_, choleskyRow_, &idum,
     620    &nrhs, integerParameters_, &msglvl, x, region, &error);
     621  if (error != 0) {
     622    printf("ERROR during solution: %d\n", error);
     623    exit(3);
     624  }
     625  delete[] x;
    633626}
    634627/* -------------------------------------------------------------------- */
     
    636629/* ..                    in pardiso solver                              */
    637630/* -------------------------------------------------------------------- */
    638     // up rows dropped
     631// up rows dropped
    639632extern "C" {
    640 int mkl_pardiso_pivot( const double*aii, double*bii, const double*eps )
     633int mkl_pardiso_pivot(const double *aii, double *bii, const double *eps)
    641634{
    642635  //if ( (*bii > *eps) || ( *bii < -*eps ) )
    643   if (*bii>*eps)
    644         return 0;
    645     if ( *bii > 0 )
    646         *bii = *eps;
    647     else
    648         *bii = -*eps;
    649     *bii = 1.0e100;
    650     return 1;
     636  if (*bii > *eps)
     637    return 0;
     638  if (*bii > 0)
     639    *bii = *eps;
     640  else
     641    *bii = -*eps;
     642  *bii = 1.0e100;
     643  return 1;
    651644}
    652645}
    653646#endif
     647
     648/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     649*/
Note: See TracChangeset for help on using the changeset viewer.