source: trunk/Clp/src/ClpCholeskyUfl.cpp @ 1197

Last change on this file since 1197 was 1197, checked in by forrest, 12 years ago

many changes to try and improve performance

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.4 KB
Line 
1#ifdef UFL_BARRIER
2// Copyright (C) 2004, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5
6
7#include "CoinPragma.hpp"
8#include "ClpCholeskyUfl.hpp"
9#include "ClpMessage.hpp"
10#include "ClpInterior.hpp"
11#include "CoinHelperFunctions.hpp"
12#include "ClpHelperFunctions.hpp"
13//#############################################################################
14// Constructors / Destructor / Assignment
15//#############################################################################
16
17//-------------------------------------------------------------------
18// Default Constructor
19//-------------------------------------------------------------------
20ClpCholeskyUfl::ClpCholeskyUfl (int denseThreshold) 
21  : ClpCholeskyBase(denseThreshold)
22{
23  type_=14;
24#ifdef CLP_USE_CHOLMOD
25  L_= NULL;
26  cholmod_start (&c_) ;
27  // Can't use supernodal as may not be positive definite
28  c_.supernodal=0; 
29#endif
30}
31
32//-------------------------------------------------------------------
33// Copy constructor
34//-------------------------------------------------------------------
35ClpCholeskyUfl::ClpCholeskyUfl (const ClpCholeskyUfl & rhs) 
36: ClpCholeskyBase(rhs)
37{
38  abort();
39}
40
41
42//-------------------------------------------------------------------
43// Destructor
44//-------------------------------------------------------------------
45ClpCholeskyUfl::~ClpCholeskyUfl ()
46{
47#ifdef CLP_USE_CHOLMOD
48  cholmod_free_factor (&L_, &c_) ;
49  cholmod_finish (&c_) ;
50#endif
51}
52
53//----------------------------------------------------------------
54// Assignment operator
55//-------------------------------------------------------------------
56ClpCholeskyUfl &
57ClpCholeskyUfl::operator=(const ClpCholeskyUfl& rhs)
58{
59  if (this != &rhs) {
60    ClpCholeskyBase::operator=(rhs);
61    abort();
62  }
63  return *this;
64}
65//-------------------------------------------------------------------
66// Clone
67//-------------------------------------------------------------------
68ClpCholeskyBase * ClpCholeskyUfl::clone() const
69{
70  return new ClpCholeskyUfl(*this);
71}
72#ifndef CLP_USE_CHOLMOD
73/* Orders rows and saves pointer to matrix.and model */
74int 
75ClpCholeskyUfl::order(ClpInterior * model) 
76{
77  int iRow;
78  model_=model;
79  if (preOrder(false,true,doKKT_))
80    return -1;
81  permuteInverse_ = new int [numberRows_];
82  permute_ = new int[numberRows_];
83  double Control[AMD_CONTROL];
84  double Info[AMD_INFO];
85
86  amd_defaults(Control);
87  //amd_control(Control);
88
89  int returnCode = amd_order(numberRows_,choleskyStart_,choleskyRow_,
90                             permute_,Control, Info);
91  delete [] choleskyRow_;
92  choleskyRow_=NULL;
93  delete [] choleskyStart_;
94  choleskyStart_=NULL;
95  //amd_info(Info);
96
97  if (returnCode!=AMD_OK) {
98    std::cout<<"AMD ordering failed"<<std::endl;
99    return 1;
100  }
101  for (iRow=0;iRow<numberRows_;iRow++) {
102    permuteInverse_[permute_[iRow]]=iRow;
103  }
104  return 0;
105}
106#else
107/* Orders rows and saves pointer to matrix.and model */
108int 
109ClpCholeskyUfl::order(ClpInterior * model) 
110{
111  numberRows_ = model->numberRows();
112  if (doKKT_) {
113    numberRows_ += numberRows_ + model->numberColumns();
114    printf("finish coding UFL KKT!\n");
115    abort();
116  }
117  rowsDropped_ = new char [numberRows_];
118  memset(rowsDropped_,0,numberRows_);
119  numberRowsDropped_=0;
120  model_=model;
121  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
122  // Space for starts
123  choleskyStart_ = new CoinBigIndex[numberRows_+1];
124  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
125  const int * columnLength = model_->clpMatrix()->getVectorLengths();
126  const int * row = model_->clpMatrix()->getIndices();
127  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
128  const int * rowLength = rowCopy_->getVectorLengths();
129  const int * column = rowCopy_->getIndices();
130  // We need two arrays for counts
131  int * which = new int [numberRows_];
132  int * used = new int[numberRows_+1];
133  CoinZeroN(used,numberRows_);
134  int iRow;
135  sizeFactor_=0;
136  for (iRow=0;iRow<numberRows_;iRow++) {
137    int number=1;
138    // make sure diagonal exists
139    which[0]=iRow;
140    used[iRow]=1;
141    if (!rowsDropped_[iRow]) {
142      CoinBigIndex startRow=rowStart[iRow];
143      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
144      for (CoinBigIndex k=startRow;k<endRow;k++) {
145        int iColumn=column[k];
146        CoinBigIndex start=columnStart[iColumn];
147        CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
148        for (CoinBigIndex j=start;j<end;j++) {
149          int jRow=row[j];
150          if (jRow>=iRow&&!rowsDropped_[jRow]) {
151            if (!used[jRow]) {
152              used[jRow]=1;
153              which[number++]=jRow;
154            }
155          }
156        }
157      }
158      sizeFactor_ += number;
159      int j;
160      for (j=0;j<number;j++)
161        used[which[j]]=0;
162    }
163  }
164  delete [] which;
165  // Now we have size - create arrays and fill in
166  try { 
167    choleskyRow_ = new int [sizeFactor_];
168  }
169  catch (...) {
170    // no memory
171    delete [] choleskyStart_;
172    choleskyStart_=NULL;
173    return -1;
174  }
175  try {
176    sparseFactor_ = new double[sizeFactor_];
177  }
178  catch (...) {
179    // no memory
180    delete [] choleskyRow_;
181    choleskyRow_=NULL;
182    delete [] choleskyStart_;
183    choleskyStart_=NULL;
184    return -1;
185  }
186 
187  sizeFactor_=0;
188  which = choleskyRow_;
189  for (iRow=0;iRow<numberRows_;iRow++) {
190    int number=1;
191    // make sure diagonal exists
192    which[0]=iRow;
193    used[iRow]=1;
194    choleskyStart_[iRow]=sizeFactor_;
195    if (!rowsDropped_[iRow]) {
196      CoinBigIndex startRow=rowStart[iRow];
197      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
198      for (CoinBigIndex k=startRow;k<endRow;k++) {
199        int iColumn=column[k];
200        CoinBigIndex start=columnStart[iColumn];
201        CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
202        for (CoinBigIndex j=start;j<end;j++) {
203          int jRow=row[j];
204          if (jRow>=iRow&&!rowsDropped_[jRow]) {
205            if (!used[jRow]) {
206              used[jRow]=1;
207              which[number++]=jRow;
208            }
209          }
210        }
211      }
212      sizeFactor_ += number;
213      int j;
214      for (j=0;j<number;j++)
215        used[which[j]]=0;
216      // Sort
217      std::sort(which,which+number);
218      // move which on
219      which += number;
220    }
221  }
222  choleskyStart_[numberRows_]=sizeFactor_;
223  delete [] used;
224  permuteInverse_ = new int [numberRows_];
225  permute_ = new int[numberRows_];
226  cholmod_sparse A ;
227  A.nrow=numberRows_;
228  A.ncol=numberRows_;
229  A.nzmax=choleskyStart_[numberRows_];
230  A.p = choleskyStart_;
231  A.i = choleskyRow_;
232  A.x=NULL;
233  A.stype=-1;
234  A.itype=CHOLMOD_INT;
235  A.xtype=CHOLMOD_PATTERN;
236  A.dtype=CHOLMOD_DOUBLE;
237  A.sorted=1;
238  A.packed=1;
239  c_.nmethods=9;
240  c_.postorder=true;
241  //c_.dbound=1.0e-20;
242  L_ = cholmod_analyze (&A, &c_) ;
243  if (c_.status) {
244    std::cout<<"CHOLMOD ordering failed"<<std::endl;
245    return 1;
246  } else {
247    printf("%g nonzeros, flop count %g\n",c_.lnz,c_.fl);
248  }
249  for (iRow=0;iRow<numberRows_;iRow++) {
250    permuteInverse_[iRow]=iRow;
251    permute_[iRow]=iRow;
252  }
253  return 0;
254}
255/* Does Symbolic factorization given permutation.
256   This is called immediately after order.  If user provides this then
257   user must provide factorize and solve.  Otherwise the default factorization is used
258   returns non-zero if not enough memory */
259int 
260ClpCholeskyUfl::symbolic()
261{
262  return 0;
263}
264/* Factorize - filling in rowsDropped and returning number dropped */
265int 
266ClpCholeskyUfl::factorize(const double * diagonal, int * rowsDropped) 
267{
268  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
269  const int * columnLength = model_->clpMatrix()->getVectorLengths();
270  const int * row = model_->clpMatrix()->getIndices();
271  const double * element = model_->clpMatrix()->getElements();
272  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
273  const int * rowLength = rowCopy_->getVectorLengths();
274  const int * column = rowCopy_->getIndices();
275  const double * elementByRow = rowCopy_->getElements();
276  int numberColumns=model_->clpMatrix()->getNumCols();
277  int iRow;
278  double * work = new double[numberRows_];
279  CoinZeroN(work,numberRows_);
280  const double * diagonalSlack = diagonal + numberColumns;
281  int newDropped=0;
282  double largest;
283  //double smallest;
284  //perturbation
285  double perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
286  perturbation=0.0;
287  perturbation=perturbation*perturbation;
288  if (perturbation>1.0) {
289#ifdef COIN_DEVELOP
290    //if (model_->model()->logLevel()&4)
291      std::cout <<"large perturbation "<<perturbation<<std::endl;
292#endif
293    perturbation=sqrt(perturbation);;
294    perturbation=1.0;
295  }
296  double delta2 = model_->delta(); // add delta*delta to diagonal
297  delta2 *= delta2;
298  for (iRow=0;iRow<numberRows_;iRow++) {
299    double * put = sparseFactor_+choleskyStart_[iRow];
300    int * which = choleskyRow_+choleskyStart_[iRow];
301    int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
302    if (!rowLength[iRow])
303      rowsDropped_[iRow]=1;
304    if (!rowsDropped_[iRow]) {
305      CoinBigIndex startRow=rowStart[iRow];
306      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
307      work[iRow] = diagonalSlack[iRow]+delta2;
308      for (CoinBigIndex k=startRow;k<endRow;k++) {
309        int iColumn=column[k];
310        if (!whichDense_||!whichDense_[iColumn]) {
311          CoinBigIndex start=columnStart[iColumn];
312          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
313          double multiplier = diagonal[iColumn]*elementByRow[k];
314          for (CoinBigIndex j=start;j<end;j++) {
315            int jRow=row[j];
316            if (jRow>=iRow&&!rowsDropped_[jRow]) {
317              double value=element[j]*multiplier;
318              work[jRow] += value;
319            }
320          }
321        }
322      }
323      int j;
324      for (j=0;j<number;j++) {
325        int jRow = which[j];
326        put[j]=work[jRow];
327        work[jRow]=0.0;
328      }
329    } else {
330      // dropped
331      int j;
332      for (j=1;j<number;j++) {
333        put[j]=0.0;
334      }
335      put[0]=1.0;
336    }
337  }
338  //check sizes
339  double largest2=maximumAbsElement(sparseFactor_,sizeFactor_);
340  largest2*=1.0e-20;
341  largest = CoinMin(largest2,1.0e-11);
342  int numberDroppedBefore=0;
343  for (iRow=0;iRow<numberRows_;iRow++) {
344    int dropped=rowsDropped_[iRow];
345    // Move to int array
346    rowsDropped[iRow]=dropped;
347    if (!dropped) {
348      CoinBigIndex start = choleskyStart_[iRow];
349      double diagonal = sparseFactor_[start];
350      if (diagonal>largest2) {
351        sparseFactor_[start]=CoinMax(diagonal,1.0e-10);
352      } else {
353        sparseFactor_[start]=CoinMax(diagonal,1.0e-10);
354        rowsDropped[iRow]=2;
355        numberDroppedBefore++;
356      } 
357    } 
358  } 
359  delete [] work;
360  cholmod_sparse A ;
361  A.nrow=numberRows_;
362  A.ncol=numberRows_;
363  A.nzmax=choleskyStart_[numberRows_];
364  A.p = choleskyStart_;
365  A.i = choleskyRow_;
366  A.x=sparseFactor_;
367  A.stype=-1;
368  A.itype=CHOLMOD_INT;
369  A.xtype=CHOLMOD_REAL;
370  A.dtype=CHOLMOD_DOUBLE;
371  A.sorted=1;
372  A.packed=1;
373  cholmod_factorize (&A, L_, &c_) ;                 /* factorize */
374  choleskyCondition_=1.0;
375  bool cleanCholesky;
376  if (model_->numberIterations()<2000) 
377    cleanCholesky=true;
378  else 
379    cleanCholesky=false;
380  if (cleanCholesky) {
381    //drop fresh makes some formADAT easier
382    //int oldDropped=numberRowsDropped_;
383    if (newDropped||numberRowsDropped_) {
384      //std::cout <<"Rank "<<numberRows_-newDropped<<" ( "<<
385      //  newDropped<<" dropped)";
386      //if (newDropped>oldDropped)
387      //std::cout<<" ( "<<newDropped-oldDropped<<" dropped this time)";
388      //std::cout<<std::endl;
389      newDropped=0;
390      for (int i=0;i<numberRows_;i++) {
391        char dropped = rowsDropped[i];
392        rowsDropped_[i]=dropped;
393        if (dropped==2) {
394          //dropped this time
395          rowsDropped[newDropped++]=i;
396          rowsDropped_[i]=0;
397        } 
398      } 
399      numberRowsDropped_=newDropped;
400      newDropped=-(2+newDropped);
401    } 
402  } else {
403    if (newDropped) {
404      newDropped=0;
405      for (int i=0;i<numberRows_;i++) {
406        char dropped = rowsDropped[i];
407        rowsDropped_[i]=dropped;
408        if (dropped==2) {
409          //dropped this time
410          rowsDropped[newDropped++]=i;
411          rowsDropped_[i]=1;
412        } 
413      } 
414    } 
415    numberRowsDropped_+=newDropped;
416    if (numberRowsDropped_&&0) {
417      std::cout <<"Rank "<<numberRows_-numberRowsDropped_<<" ( "<<
418          numberRowsDropped_<<" dropped)";
419      if (newDropped) {
420        std::cout<<" ( "<<newDropped<<" dropped this time)";
421      } 
422      std::cout<<std::endl;
423    } 
424  }
425  status_=0;
426  return newDropped;
427}
428/* Uses factorization to solve. */
429void 
430ClpCholeskyUfl::solve (double * region) 
431{
432  cholmod_dense *x, *b;
433  b = cholmod_allocate_dense (numberRows_, 1, numberRows_, CHOLMOD_REAL, &c_) ; 
434  CoinMemcpyN(region,numberRows_,(double *) b->x);
435  x = cholmod_solve (CHOLMOD_A, L_, b, &c_) ;
436  CoinMemcpyN((double *) x->x,numberRows_,region);
437  cholmod_free_dense (&x, &c_) ;
438  cholmod_free_dense (&b, &c_) ;
439}
440#endif
441#endif
Note: See TracBrowser for help on using the repository browser.