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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

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