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

Last change on this file since 1034 was 1034, checked in by forrest, 13 years ago

moving branches/devel to trunk

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