source: trunk/Clp/src/ClpCholeskyPardiso.cpp @ 2385

Last change on this file since 2385 was 2385, checked in by unxusr, 4 months ago

formatting

File size: 24.7 KB
Line 
1/* $Id: ClpCholeskyPardiso.cpp 1723 2011-04-17 15:07:10Z forrest $ */
2// Copyright (C) 2003, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5#ifdef PARDISO_BARRIER
6
7#include "CoinPragma.hpp"
8#include "CoinHelperFunctions.hpp"
9#include "ClpHelperFunctions.hpp"
10
11#include "ClpInterior.hpp"
12#include "ClpCholeskyPardiso.hpp"
13#include "ClpCholeskyDense.hpp"
14#include "ClpMessage.hpp"
15
16//#############################################################################
17// Constructors / Destructor / Assignment
18//#############################################################################
19
20//-------------------------------------------------------------------
21// Default Constructor
22//-------------------------------------------------------------------
23ClpCholeskyPardiso::ClpCholeskyPardiso(int denseThreshold)
24  : ClpCholeskyBase(denseThreshold)
25{
26  type_ = 17;
27}
28
29//-------------------------------------------------------------------
30// Copy constructor
31//-------------------------------------------------------------------
32ClpCholeskyPardiso::ClpCholeskyPardiso(const ClpCholeskyPardiso &rhs)
33  : ClpCholeskyBase(rhs)
34{
35}
36
37//-------------------------------------------------------------------
38// Destructor
39//-------------------------------------------------------------------
40ClpCholeskyPardiso::~ClpCholeskyPardiso()
41{
42  /* -------------------------------------------------------------------- */
43  /* .. Termination and release of memory. */
44  /* -------------------------------------------------------------------- */
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;
51  /* Auxiliary variables. */
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);
58}
59
60//----------------------------------------------------------------
61// Assignment operator
62//-------------------------------------------------------------------
63ClpCholeskyPardiso &
64ClpCholeskyPardiso::operator=(const ClpCholeskyPardiso &rhs)
65{
66  if (this != &rhs) {
67    ClpCholeskyBase::operator=(rhs);
68  }
69  return *this;
70}
71//-------------------------------------------------------------------
72// Clone
73//-------------------------------------------------------------------
74ClpCholeskyBase *ClpCholeskyPardiso::clone() const
75{
76  return new ClpCholeskyPardiso(*this);
77}
78/* Orders rows and saves pointer to matrix.and model */
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            }
169          }
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            }
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;
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;
290#define ALWAYS_ORDER
291#ifdef ALWAYS_ORDER
292  error = 0;
293#else
294  phase = 11;
295  PARDISO(voidParameters, &maxfct, &mnum, &mtype, &phase,
296    &numberRows_, sparseFactor_,
297    choleskyStart_, choleskyRow_, &idum,
298    &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error);
299#endif
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;
313}
314/* Does Symbolic factorization given permutation.
315   This is called immediately after order.  If user provides this then
316   user must provide factorize and solve.  Otherwise the default factorization is used
317   returns non-zero if not enough memory */
318int ClpCholeskyPardiso::symbolic()
319{
320  return 0;
321}
322/* Factorize - filling in rowsDropped and returning number dropped */
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) {
346#ifdef COIN_DEVELOP
347    //if (model_->model()->logLevel()&4)
348    std::cout << "large perturbation " << perturbation << std::endl;
349#endif
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            }
406          }
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];
551          }
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;
600}
601/* Uses factorization to solve. */
602void ClpCholeskyPardiso::solve(double *region)
603{
604  /* -------------------------------------------------------------------- */
605  /* .. Back substitution and iterative refinement. */
606  /* -------------------------------------------------------------------- */
607  int phase = 33;
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;
615  /* Auxiliary variables. */
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;
626}
627/* -------------------------------------------------------------------- */
628/* .. mkl_pardiso_pivot: This function replace appropriate function     */
629/* ..                    in pardiso solver                              */
630/* -------------------------------------------------------------------- */
631// up rows dropped
632extern "C" {
633int mkl_pardiso_pivot(const double *aii, double *bii, const double *eps)
634{
635  //if ( (*bii > *eps) || ( *bii < -*eps ) )
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;
644}
645}
646#endif
647
648/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
649*/
Note: See TracBrowser for help on using the repository browser.