source: trunk/Clp/src/ClpCholeskyWssmpKKT.cpp @ 2470

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

formatting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 18.9 KB
Line 
1/* $Id: ClpCholeskyWssmpKKT.cpp 2385 2019-01-06 19:43:06Z stefan $ */
2// Copyright (C) 2004, 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
6#include "CoinPragma.hpp"
7#include "CoinHelperFunctions.hpp"
8#include "ClpHelperFunctions.hpp"
9
10#include "ClpInterior.hpp"
11#include "ClpCholeskyWssmpKKT.hpp"
12#include "ClpQuadraticObjective.hpp"
13#include "ClpMessage.hpp"
14
15//#############################################################################
16// Constructors / Destructor / Assignment
17//#############################################################################
18
19//-------------------------------------------------------------------
20// Default Constructor
21//-------------------------------------------------------------------
22ClpCholeskyWssmpKKT::ClpCholeskyWssmpKKT(int denseThreshold)
23  : ClpCholeskyBase(denseThreshold)
24{
25  type_ = 21;
26}
27
28//-------------------------------------------------------------------
29// Copy constructor
30//-------------------------------------------------------------------
31ClpCholeskyWssmpKKT::ClpCholeskyWssmpKKT(const ClpCholeskyWssmpKKT &rhs)
32  : ClpCholeskyBase(rhs)
33{
34}
35
36//-------------------------------------------------------------------
37// Destructor
38//-------------------------------------------------------------------
39ClpCholeskyWssmpKKT::~ClpCholeskyWssmpKKT()
40{
41}
42
43//----------------------------------------------------------------
44// Assignment operator
45//-------------------------------------------------------------------
46ClpCholeskyWssmpKKT &
47ClpCholeskyWssmpKKT::operator=(const ClpCholeskyWssmpKKT &rhs)
48{
49  if (this != &rhs) {
50    ClpCholeskyBase::operator=(rhs);
51  }
52  return *this;
53}
54//-------------------------------------------------------------------
55// Clone
56//-------------------------------------------------------------------
57ClpCholeskyBase *ClpCholeskyWssmpKKT::clone() const
58{
59  return new ClpCholeskyWssmpKKT(*this);
60}
61// At present I can't get wssmp to work as my libraries seem to be out of sync
62// so I have linked in ekkwssmp which is an older version
63#ifndef USE_EKKWSSMP
64extern "C" {
65void F77_FUNC(wsetmaxthrds, WSETMAXTHRDS)(const int *NTHREADS);
66
67void F77_FUNC(wssmp, WSSMP)(const int *N, const int *IA,
68  const int *JA, const double *AVALS,
69  double *DIAG, int *PERM,
70  int *INVP, double *B,
71  const int *LDB, const int *NRHS,
72  double *AUX, const int *NAUX,
73  int *MRP, int *IPARM,
74  double *DPARM);
75void F77_FUNC_(wsmp_clear, WSMP_CLEAR)(void);
76}
77#else
78/* minimum needed for user */
79typedef struct EKKModel EKKModel;
80typedef struct EKKContext EKKContext;
81
82extern "C" {
83EKKContext *ekk_initializeContext();
84void ekk_endContext(EKKContext *context);
85EKKModel *ekk_newModel(EKKContext *env, const char *name);
86int ekk_deleteModel(EKKModel *model);
87}
88static EKKModel *model = NULL;
89static EKKContext *context = NULL;
90extern "C" void ekkwssmp(EKKModel *, int *n,
91  int *columnStart, int *rowIndex, double *element,
92  double *diagonal, int *perm, int *invp,
93  double *rhs, int *ldb, int *nrhs,
94  double *aux, int *naux,
95  int *mrp, int *iparm, double *dparm);
96static void F77_FUNC(wssmp, WSSMP)(int *n, int *ia, int *ja,
97  double *avals, double *diag, int *perm, int *invp,
98  double *b, int *ldb, int *nrhs, double *aux, int *naux, int *mrp, int *iparm, double *dparm)
99{
100  if (!context) {
101    /* initialize OSL environment */
102    context = ekk_initializeContext();
103    model = ekk_newModel(context, "");
104  }
105  ekkwssmp(model, n, ia, ja,
106    avals, diag, perm, invp,
107    b, ldb, nrhs, aux,
108    naux, mrp, iparm, dparm);
109  //ekk_deleteModel(model);
110  //ekk_endContext(context);
111}
112#endif
113
114/* Orders rows and saves pointer to model */
115int ClpCholeskyWssmpKKT::order(ClpInterior *model)
116{
117  int numberRowsModel = model->numberRows();
118  int numberColumns = model->numberColumns();
119  int numberTotal = numberColumns + numberRowsModel;
120  numberRows_ = 2 * numberRowsModel + numberColumns;
121  rowsDropped_ = new char[numberRows_];
122  memset(rowsDropped_, 0, numberRows_);
123  numberRowsDropped_ = 0;
124  model_ = model;
125  CoinPackedMatrix *quadratic = NULL;
126  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(model_->objectiveAsObject()));
127  if (quadraticObj)
128    quadratic = quadraticObj->quadraticObjective();
129  int numberElements = model_->clpMatrix()->getNumElements();
130  numberElements = numberElements + 2 * numberRowsModel + numberTotal;
131  if (quadratic)
132    numberElements += quadratic->getNumElements();
133  // Space for starts
134  choleskyStart_ = new CoinBigIndex[numberRows_ + 1];
135  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
136  const int *columnLength = model_->clpMatrix()->getVectorLengths();
137  const int *row = model_->clpMatrix()->getIndices();
138  //const double * element = model_->clpMatrix()->getElements();
139  // Now we have size - create arrays and fill in
140  try {
141    choleskyRow_ = new int[numberElements];
142  } catch (...) {
143    // no memory
144    delete[] choleskyStart_;
145    choleskyStart_ = NULL;
146    return -1;
147  }
148  try {
149    sparseFactor_ = new double[numberElements];
150  } catch (...) {
151    // no memory
152    delete[] choleskyRow_;
153    choleskyRow_ = NULL;
154    delete[] choleskyStart_;
155    choleskyStart_ = NULL;
156    return -1;
157  }
158  int iRow, iColumn;
159
160  sizeFactor_ = 0;
161  // matrix
162  if (!quadratic) {
163    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
164      choleskyStart_[iColumn] = sizeFactor_;
165      choleskyRow_[sizeFactor_++] = iColumn;
166      CoinBigIndex start = columnStart[iColumn];
167      CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
168      for (CoinBigIndex j = start; j < end; j++) {
169        choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
170      }
171    }
172  } else {
173    // Quadratic
174    const int *columnQuadratic = quadratic->getIndices();
175    const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
176    const int *columnQuadraticLength = quadratic->getVectorLengths();
177    //const double * quadraticElement = quadratic->getElements();
178    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
179      choleskyStart_[iColumn] = sizeFactor_;
180      choleskyRow_[sizeFactor_++] = iColumn;
181      for (CoinBigIndex j = columnQuadraticStart[iColumn];
182           j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
183        int jColumn = columnQuadratic[j];
184        if (jColumn > iColumn)
185          choleskyRow_[sizeFactor_++] = jColumn;
186      }
187      CoinBigIndex start = columnStart[iColumn];
188      CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
189      for (CoinBigIndex j = start; j < end; j++) {
190        choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
191      }
192    }
193  }
194  // slacks
195  for (; iColumn < numberTotal; iColumn++) {
196    choleskyStart_[iColumn] = sizeFactor_;
197    choleskyRow_[sizeFactor_++] = iColumn;
198    choleskyRow_[sizeFactor_++] = iColumn - numberColumns + numberTotal;
199  }
200  // Transpose - nonzero diagonal (may regularize)
201  for (iRow = 0; iRow < numberRowsModel; iRow++) {
202    choleskyStart_[iRow + numberTotal] = sizeFactor_;
203    // diagonal
204    choleskyRow_[sizeFactor_++] = iRow + numberTotal;
205  }
206  choleskyStart_[numberRows_] = sizeFactor_;
207  permuteInverse_ = new int[numberRows_];
208  permute_ = new int[numberRows_];
209  integerParameters_[0] = 0;
210  int i0 = 0;
211  int i1 = 1;
212#ifndef USE_EKKWSSMP
213  int i2 = 1;
214  if (model->numberThreads() <= 0)
215    i2 = 1;
216  else
217    i2 = model->numberThreads();
218  F77_FUNC(wsetmaxthrds, WSETMAXTHRDS)
219  (&i2);
220#endif
221  F77_FUNC(wssmp, WSSMP)
222  (&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
223    NULL, permute_, permuteInverse_, 0, &numberRows_, &i1,
224    NULL, &i0, NULL, integerParameters_, doubleParameters_);
225  integerParameters_[1] = 1; //order and symbolic
226  integerParameters_[2] = 2;
227  integerParameters_[3] = 0; //CSR
228  integerParameters_[4] = 0; //C style
229  integerParameters_[13] = 1; //reuse initial factorization space
230  integerParameters_[15 + 0] = 1; //ordering
231  integerParameters_[15 + 1] = 0;
232  integerParameters_[15 + 2] = 1;
233  integerParameters_[15 + 3] = 0;
234  integerParameters_[15 + 4] = 1;
235  doubleParameters_[10] = 1.0e-20;
236  doubleParameters_[11] = 1.0e-15;
237#if 1
238  integerParameters_[1] = 2; //just symbolic
239  for (int iRow = 0; iRow < numberRows_; iRow++) {
240    permuteInverse_[iRow] = iRow;
241    permute_[iRow] = iRow;
242  }
243#endif
244  F77_FUNC(wssmp, WSSMP)
245  (&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
246    NULL, permute_, permuteInverse_, NULL, &numberRows_, &i1,
247    NULL, &i0, NULL, integerParameters_, doubleParameters_);
248  //std::cout<<"Ordering and symbolic factorization took "<<doubleParameters_[0]<<std::endl;
249  if (integerParameters_[63]) {
250    std::cout << "wssmp returning error code of " << integerParameters_[63] << std::endl;
251    return 1;
252  }
253  std::cout << integerParameters_[23] << " elements in sparse Cholesky" << std::endl;
254  if (!integerParameters_[23]) {
255    for (int iRow = 0; iRow < numberRows_; iRow++) {
256      permuteInverse_[iRow] = iRow;
257      permute_[iRow] = iRow;
258    }
259    std::cout << "wssmp says no elements - fully dense? - switching to dense" << std::endl;
260    integerParameters_[1] = 2;
261    integerParameters_[2] = 2;
262    integerParameters_[7] = 1; // no permute
263    F77_FUNC(wssmp, WSSMP)
264    (&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
265      NULL, permute_, permuteInverse_, NULL, &numberRows_, &i1,
266      NULL, &i0, NULL, integerParameters_, doubleParameters_);
267    std::cout << integerParameters_[23] << " elements in dense Cholesky" << std::endl;
268  }
269  return 0;
270}
271/* Does Symbolic factorization given permutation.
272   This is called immediately after order.  If user provides this then
273   user must provide factorize and solve.  Otherwise the default factorization is used
274   returns non-zero if not enough memory */
275int ClpCholeskyWssmpKKT::symbolic()
276{
277  return 0;
278}
279/* Factorize - filling in rowsDropped and returning number dropped */
280int ClpCholeskyWssmpKKT::factorize(const double *diagonal, int *rowsDropped)
281{
282  int numberRowsModel = model_->numberRows();
283  int numberColumns = model_->numberColumns();
284  int numberTotal = numberColumns + numberRowsModel;
285  int newDropped = 0;
286  double largest = 0.0;
287  double smallest;
288  //perturbation
289  double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
290  perturbation = perturbation * perturbation;
291  if (perturbation > 1.0) {
292#ifdef COIN_DEVELOP
293    //if (model_->model()->logLevel()&4)
294    std::cout << "large perturbation " << perturbation << std::endl;
295#endif
296    perturbation = sqrt(perturbation);
297    ;
298    perturbation = 1.0;
299  }
300  // need to recreate every time
301  int iRow, iColumn;
302  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
303  const int *columnLength = model_->clpMatrix()->getVectorLengths();
304  const int *row = model_->clpMatrix()->getIndices();
305  const double *element = model_->clpMatrix()->getElements();
306
307  CoinBigIndex numberElements = 0;
308  CoinPackedMatrix *quadratic = NULL;
309  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(model_->objectiveAsObject()));
310  if (quadraticObj)
311    quadratic = quadraticObj->quadraticObjective();
312  // matrix
313  if (!quadratic) {
314    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
315      choleskyStart_[iColumn] = numberElements;
316      double value = diagonal[iColumn];
317      if (fabs(value) > 1.0e-100) {
318        value = 1.0 / value;
319        largest = CoinMax(largest, fabs(value));
320        sparseFactor_[numberElements] = -value;
321        choleskyRow_[numberElements++] = iColumn;
322        CoinBigIndex start = columnStart[iColumn];
323        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
324        for (CoinBigIndex j = start; j < end; j++) {
325          choleskyRow_[numberElements] = row[j] + numberTotal;
326          sparseFactor_[numberElements++] = element[j];
327          largest = CoinMax(largest, fabs(element[j]));
328        }
329      } else {
330        sparseFactor_[numberElements] = -1.0e100;
331        choleskyRow_[numberElements++] = iColumn;
332      }
333    }
334  } else {
335    // Quadratic
336    const int *columnQuadratic = quadratic->getIndices();
337    const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
338    const int *columnQuadraticLength = quadratic->getVectorLengths();
339    const double *quadraticElement = quadratic->getElements();
340    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
341      choleskyStart_[iColumn] = numberElements;
342      CoinBigIndex savePosition = numberElements;
343      choleskyRow_[numberElements++] = iColumn;
344      double value = diagonal[iColumn];
345      if (fabs(value) > 1.0e-100) {
346        value = 1.0 / value;
347        for (CoinBigIndex j = columnQuadraticStart[iColumn];
348             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
349          int jColumn = columnQuadratic[j];
350          if (jColumn > iColumn) {
351            sparseFactor_[numberElements] = -quadraticElement[j];
352            choleskyRow_[numberElements++] = jColumn;
353          } else if (iColumn == jColumn) {
354            value += quadraticElement[j];
355          }
356        }
357        largest = CoinMax(largest, fabs(value));
358        sparseFactor_[savePosition] = -value;
359        CoinBigIndex start = columnStart[iColumn];
360        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
361        for (CoinBigIndex j = start; j < end; j++) {
362          choleskyRow_[numberElements] = row[j] + numberTotal;
363          sparseFactor_[numberElements++] = element[j];
364          largest = CoinMax(largest, fabs(element[j]));
365        }
366      } else {
367        value = 1.0e100;
368        sparseFactor_[savePosition] = -value;
369      }
370    }
371  }
372  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
373    assert(sparseFactor_[choleskyStart_[iColumn]] < 0.0);
374  }
375  // slacks
376  for (iColumn = numberColumns; iColumn < numberTotal; iColumn++) {
377    choleskyStart_[iColumn] = numberElements;
378    double value = diagonal[iColumn];
379    if (fabs(value) > 1.0e-100) {
380      value = 1.0 / value;
381      largest = CoinMax(largest, fabs(value));
382    } else {
383      value = 1.0e100;
384    }
385    sparseFactor_[numberElements] = -value;
386    choleskyRow_[numberElements++] = iColumn;
387    choleskyRow_[numberElements] = iColumn - numberColumns + numberTotal;
388    sparseFactor_[numberElements++] = -1.0;
389  }
390  // Finish diagonal
391  double delta2 = model_->delta(); // add delta*delta to bottom
392  delta2 *= delta2;
393  for (iRow = 0; iRow < numberRowsModel; iRow++) {
394    choleskyStart_[iRow + numberTotal] = numberElements;
395    choleskyRow_[numberElements] = iRow + numberTotal;
396    sparseFactor_[numberElements++] = delta2;
397  }
398  choleskyStart_[numberRows_] = numberElements;
399  int i1 = 1;
400  int i0 = 0;
401  integerParameters_[1] = 3;
402  integerParameters_[2] = 3;
403  integerParameters_[10] = 2;
404  //integerParameters_[11]=1;
405  integerParameters_[12] = 2;
406  // LDLT
407  integerParameters_[30] = 1;
408  doubleParameters_[20] = 1.0e100;
409  double largest2 = largest * 1.0e-20;
410  largest = CoinMin(largest2, 1.0e-11);
411  doubleParameters_[10] = CoinMax(1.0e-20, largest);
412  if (doubleParameters_[10] > 1.0e-3)
413    integerParameters_[9] = 1;
414  else
415    integerParameters_[9] = 0;
416#ifndef WSMP
417  // Set up LDL cutoff
418  integerParameters_[34] = numberTotal;
419  doubleParameters_[20] = 1.0e-15;
420  doubleParameters_[34] = 1.0e-12;
421  //printf("tol is %g\n",doubleParameters_[10]);
422  //doubleParameters_[10]=1.0e-17;
423#endif
424  int *rowsDropped2 = new int[numberRows_];
425  CoinZeroN(rowsDropped2, numberRows_);
426  F77_FUNC(wssmp, WSSMP)
427  (&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
428    NULL, permute_, permuteInverse_, NULL, &numberRows_, &i1,
429    NULL, &i0, rowsDropped2, integerParameters_, doubleParameters_);
430  //std::cout<<"factorization took "<<doubleParameters_[0]<<std::endl;
431  if (integerParameters_[9]) {
432    std::cout << "scaling applied" << std::endl;
433  }
434  newDropped = integerParameters_[20];
435#if 1
436  // Should save adjustments in ..R_
437  int n1 = 0, n2 = 0;
438  double *primalR = model_->primalR();
439  double *dualR = model_->dualR();
440  for (iRow = 0; iRow < numberTotal; iRow++) {
441    if (rowsDropped2[iRow]) {
442      n1++;
443      //printf("row region1 %d dropped\n",iRow);
444      //rowsDropped_[iRow]=1;
445      rowsDropped_[iRow] = 0;
446      primalR[iRow] = doubleParameters_[20];
447    } else {
448      rowsDropped_[iRow] = 0;
449      primalR[iRow] = 0.0;
450    }
451  }
452  for (; iRow < numberRows_; iRow++) {
453    if (rowsDropped2[iRow]) {
454      n2++;
455      //printf("row region2 %d dropped\n",iRow);
456      //rowsDropped_[iRow]=1;
457      rowsDropped_[iRow] = 0;
458      dualR[iRow - numberTotal] = doubleParameters_[34];
459    } else {
460      rowsDropped_[iRow] = 0;
461      dualR[iRow - numberTotal] = 0.0;
462    }
463  }
464  //printf("%d rows dropped in region1, %d in region2\n",n1,n2);
465#endif
466  delete[] rowsDropped2;
467  //if (integerParameters_[20])
468  //std::cout<<integerParameters_[20]<<" rows dropped"<<std::endl;
469  largest = doubleParameters_[3];
470  smallest = doubleParameters_[4];
471  if (model_->messageHandler()->logLevel() > 1)
472    std::cout << "Cholesky - largest " << largest << " smallest " << smallest << std::endl;
473  choleskyCondition_ = largest / smallest;
474  if (integerParameters_[63] < 0)
475    return -1; // out of memory
476  status_ = 0;
477  return 0;
478}
479/* Uses factorization to solve. */
480void ClpCholeskyWssmpKKT::solve(double *region)
481{
482  abort();
483}
484/* Uses factorization to solve. */
485void ClpCholeskyWssmpKKT::solveKKT(double *region1, double *region2, const double *diagonal,
486  double diagonalScaleFactor)
487{
488  int numberRowsModel = model_->numberRows();
489  int numberColumns = model_->numberColumns();
490  int numberTotal = numberColumns + numberRowsModel;
491  double *array = new double[numberRows_];
492  CoinMemcpyN(region1, numberTotal, array);
493  CoinMemcpyN(region2, numberRowsModel, array + numberTotal);
494  int i1 = 1;
495  int i0 = 0;
496  integerParameters_[1] = 4;
497  integerParameters_[2] = 4;
498#if 0
499     integerParameters_[5] = 3;
500     doubleParameters_[5] = 1.0e-10;
501     integerParameters_[6] = 6;
502#endif
503  F77_FUNC(wssmp, WSSMP)
504  (&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
505    NULL, permute_, permuteInverse_, array, &numberRows_, &i1,
506    NULL, &i0, NULL, integerParameters_, doubleParameters_);
507#if 0
508     int iRow;
509     for (iRow = 0; iRow < numberTotal; iRow++) {
510          if (rowsDropped_[iRow] && fabs(array[iRow]) > 1.0e-8) {
511               printf("row region1 %d dropped %g\n", iRow, array[iRow]);
512          }
513     }
514     for (; iRow < numberRows_; iRow++) {
515          if (rowsDropped_[iRow] && fabs(array[iRow]) > 1.0e-8) {
516               printf("row region2 %d dropped %g\n", iRow, array[iRow]);
517          }
518     }
519#endif
520  CoinMemcpyN(array + numberTotal, numberRowsModel, region2);
521#if 1
522  CoinMemcpyN(array, numberTotal, region1);
523#else
524  multiplyAdd(region2, numberRowsModel, -1.0, array + numberColumns, 0.0);
525  CoinZeroN(array, numberColumns);
526  model_->clpMatrix()->transposeTimes(1.0, region2, array);
527  for (int iColumn = 0; iColumn < numberTotal; iColumn++)
528    region1[iColumn] = diagonal[iColumn] * (array[iColumn] - region1[iColumn]);
529#endif
530  delete[] array;
531#if 0
532     if (integerParameters_[5]) {
533          std::cout << integerParameters_[5] << " refinements ";
534     }
535     std::cout << doubleParameters_[6] << std::endl;
536#endif
537}
538
539/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
540*/
Note: See TracBrowser for help on using the repository browser.