source: tags/move-to-subversion/ClpCholeskyWssmpKKT.cpp @ 1355

Last change on this file since 1355 was 740, checked in by forrest, 14 years ago

makes some problems faster

  • 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    //if (model_->model()->logLevel()&4)
314    std::cout <<"large perturbation "<<perturbation<<std::endl;
315    perturbation=sqrt(perturbation);;
316    perturbation=1.0;
317  }
318  // need to recreate every time
319  int iRow,iColumn;
320  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
321  const int * columnLength = model_->clpMatrix()->getVectorLengths();
322  const int * row = model_->clpMatrix()->getIndices();
323  const double * element = model_->clpMatrix()->getElements();
324 
325  CoinBigIndex numberElements=0;
326  CoinPackedMatrix * quadratic = NULL;
327  ClpQuadraticObjective * quadraticObj = 
328    (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
329  if (quadraticObj) 
330    quadratic = quadraticObj->quadraticObjective();
331  // matrix
332  if (!quadratic) {
333    for (iColumn=0;iColumn<numberColumns;iColumn++) {
334      choleskyStart_[iColumn]=numberElements;
335      double value = diagonal[iColumn];
336      if (fabs(value)>1.0e-100) {
337        value = 1.0/value;
338        largest = CoinMax(largest,fabs(value));
339        sparseFactor_[numberElements] = -value;
340        choleskyRow_[numberElements++]=iColumn;
341        CoinBigIndex start=columnStart[iColumn];
342        CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
343        for (CoinBigIndex j=start;j<end;j++) {
344          choleskyRow_[numberElements]=row[j]+numberTotal;
345          sparseFactor_[numberElements++]=element[j];
346          largest = CoinMax(largest,fabs(element[j]));
347        }
348      } else {
349        sparseFactor_[numberElements] = -1.0e100;
350        choleskyRow_[numberElements++]=iColumn;
351      }
352    }
353  } else {
354    // Quadratic
355    const int * columnQuadratic = quadratic->getIndices();
356    const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
357    const int * columnQuadraticLength = quadratic->getVectorLengths();
358    const double * quadraticElement = quadratic->getElements();
359    for (iColumn=0;iColumn<numberColumns;iColumn++) {
360      choleskyStart_[iColumn]=numberElements;
361      CoinBigIndex savePosition = numberElements;
362      choleskyRow_[numberElements++]=iColumn;
363      double value = diagonal[iColumn];
364      if (fabs(value)>1.0e-100) {
365        value = 1.0/value;
366        for (CoinBigIndex j=columnQuadraticStart[iColumn];
367             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
368          int jColumn = columnQuadratic[j];
369          if (jColumn>iColumn) {
370            sparseFactor_[numberElements]=-quadraticElement[j];
371            choleskyRow_[numberElements++]=jColumn;
372          } else if (iColumn==jColumn) {
373            value += quadraticElement[j];
374          }
375        }
376        largest = CoinMax(largest,fabs(value));
377        sparseFactor_[savePosition] = -value;
378        CoinBigIndex start=columnStart[iColumn];
379        CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
380        for (CoinBigIndex j=start;j<end;j++) {
381          choleskyRow_[numberElements]=row[j]+numberTotal;
382          sparseFactor_[numberElements++]=element[j];
383          largest = CoinMax(largest,fabs(element[j]));
384        }
385      } else {
386        value = 1.0e100;
387        sparseFactor_[savePosition] = -value;
388      }
389    }
390  }
391  for (iColumn=0;iColumn<numberColumns;iColumn++) {
392    assert (sparseFactor_[choleskyStart_[iColumn]]<0.0);
393  }
394  // slacks
395  for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) {
396    choleskyStart_[iColumn]=numberElements;
397    double value = diagonal[iColumn];
398    if (fabs(value)>1.0e-100) {
399      value = 1.0/value;
400      largest = CoinMax(largest,fabs(value));
401    } else {
402      value = 1.0e100;
403    }
404    sparseFactor_[numberElements] = -value;
405    choleskyRow_[numberElements++]=iColumn;
406    choleskyRow_[numberElements]=iColumn-numberColumns+numberTotal;
407    sparseFactor_[numberElements++]=-1.0;
408  }
409  // Finish diagonal
410  double delta2 = model_->delta(); // add delta*delta to bottom
411  delta2 *= delta2;
412  for (iRow=0;iRow<numberRowsModel;iRow++) {
413    choleskyStart_[iRow+numberTotal]=numberElements;
414    choleskyRow_[numberElements]=iRow+numberTotal;
415    sparseFactor_[numberElements++]=delta2;
416  }
417  choleskyStart_[numberRows_]=numberElements;
418  int i1=1;
419  int i0=0;
420  integerParameters_[1]=3;
421  integerParameters_[2]=3;
422  integerParameters_[10]=2;
423  //integerParameters_[11]=1;
424  integerParameters_[12]=2;
425  // LDLT
426  integerParameters_[30]=1;
427  doubleParameters_[20]=1.0e100;
428  double largest2= largest*1.0e-20;
429  largest = CoinMin(largest2,1.0e-11);
430  doubleParameters_[10]=CoinMax(1.0e-20,largest);
431  if (doubleParameters_[10]>1.0e-3)
432    integerParameters_[9]=1;
433  else
434    integerParameters_[9]=0;
435#ifndef WSMP
436  // Set up LDL cutoff
437  integerParameters_[34]=numberTotal;
438  doubleParameters_[20]=1.0e-15;
439  doubleParameters_[34]=1.0e-12;
440  //printf("tol is %g\n",doubleParameters_[10]);
441  //doubleParameters_[10]=1.0e-17;
442#endif
443  int * rowsDropped2 = new int[numberRows_];
444  CoinZeroN(rowsDropped2,numberRows_);
445  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
446        NULL,permute_,permuteInverse_,NULL,&numberRows_,&i1,
447        NULL,&i0,rowsDropped2,integerParameters_,doubleParameters_);
448   //std::cout<<"factorization took "<<doubleParameters_[0]<<std::endl;
449  if (integerParameters_[9]) {
450    std::cout<<"scaling applied"<<std::endl;
451  } 
452  newDropped=integerParameters_[20];
453#if 1
454  // Should save adjustments in ..R_
455  int n1=0,n2=0;
456  double * primalR = model_->primalR();
457  double * dualR = model_->dualR();
458  for (iRow=0;iRow<numberTotal;iRow++) { 
459    if (rowsDropped2[iRow]) {
460      n1++;
461      //printf("row region1 %d dropped\n",iRow);
462      //rowsDropped_[iRow]=1;
463      rowsDropped_[iRow]=0;
464      primalR[iRow]=doubleParameters_[20];
465    } else {
466      rowsDropped_[iRow]=0;
467      primalR[iRow]=0.0;
468    }
469  }
470  for (;iRow<numberRows_;iRow++) {
471    if (rowsDropped2[iRow]) {
472      n2++;
473      //printf("row region2 %d dropped\n",iRow);
474      //rowsDropped_[iRow]=1;
475      rowsDropped_[iRow]=0;
476      dualR[iRow-numberTotal]=doubleParameters_[34];
477    } else {
478      rowsDropped_[iRow]=0;
479      dualR[iRow-numberTotal]=0.0;
480    }
481  }
482  //printf("%d rows dropped in region1, %d in region2\n",n1,n2);
483#endif
484  delete [] rowsDropped2;
485  //if (integerParameters_[20])
486  //std::cout<<integerParameters_[20]<<" rows dropped"<<std::endl;
487  largest=doubleParameters_[3];
488  smallest=doubleParameters_[4];
489  if (model_->messageHandler()->logLevel()>1) 
490    std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
491  choleskyCondition_=largest/smallest;
492  if (integerParameters_[63]<0)
493    return -1; // out of memory
494  status_=0;
495  return 0;
496}
497/* Uses factorization to solve. */
498void 
499ClpCholeskyWssmpKKT::solve (double * region) 
500{
501  abort();
502}
503/* Uses factorization to solve. */
504void 
505ClpCholeskyWssmpKKT::solveKKT (double * region1, double * region2, const double * diagonal,
506                               double diagonalScaleFactor) 
507{
508  int numberRowsModel = model_->numberRows();
509  int numberColumns = model_->numberColumns();
510  int numberTotal = numberColumns + numberRowsModel;
511  double * array = new double [numberRows_];
512  CoinMemcpyN(region1,numberTotal,array);
513  CoinMemcpyN(region2,numberRowsModel,array+numberTotal);
514  int i1=1;
515  int i0=0;
516  integerParameters_[1]=4;
517  integerParameters_[2]=4;
518#if 0
519  integerParameters_[5]=3;
520  doubleParameters_[5]=1.0e-10;
521  integerParameters_[6]=6;
522#endif
523  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
524        NULL,permute_,permuteInverse_,array,&numberRows_,&i1,
525        NULL,&i0,NULL,integerParameters_,doubleParameters_);
526#if 1
527  int iRow;
528  for (iRow=0;iRow<numberTotal;iRow++) { 
529    if (rowsDropped_[iRow]&&fabs(array[iRow])>1.0e-8) {
530      printf("row region1 %d dropped %g\n",iRow,array[iRow]);
531    }
532  }
533  for (;iRow<numberRows_;iRow++) {
534    if (rowsDropped_[iRow]&&fabs(array[iRow])>1.0e-8) {
535      printf("row region2 %d dropped %g\n",iRow,array[iRow]);
536    }
537  }
538#endif
539  CoinMemcpyN(array+numberTotal,numberRowsModel,region2);
540#if 1
541  CoinMemcpyN(array,numberTotal,region1);
542#else
543  multiplyAdd(region2,numberRowsModel,-1.0,array+numberColumns,0.0);
544  CoinZeroN(array,numberColumns);
545  model_->clpMatrix()->transposeTimes(1.0,region2,array);
546  for (int iColumn=0;iColumn<numberTotal;iColumn++) 
547    region1[iColumn] = diagonal[iColumn]*(array[iColumn]-region1[iColumn]);
548#endif
549  delete [] array;
550#if 0
551  if (integerParameters_[5]) {
552    std::cout<<integerParameters_[5]<<" refinements ";
553  }
554  std::cout<<doubleParameters_[6]<<std::endl;
555#endif
556}
557#endif
Note: See TracBrowser for help on using the repository browser.