source: stable/1.3/Clp/src/ClpCholeskyUfl.cpp @ 979

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

switch off supernodal

  • 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    //if (model_->model()->logLevel()&4)
290      std::cout <<"large perturbation "<<perturbation<<std::endl;
291    perturbation=sqrt(perturbation);;
292    perturbation=1.0;
293  }
294  double delta2 = model_->delta(); // add delta*delta to diagonal
295  delta2 *= delta2;
296  for (iRow=0;iRow<numberRows_;iRow++) {
297    double * put = sparseFactor_+choleskyStart_[iRow];
298    int * which = choleskyRow_+choleskyStart_[iRow];
299    int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
300    if (!rowLength[iRow])
301      rowsDropped_[iRow]=1;
302    if (!rowsDropped_[iRow]) {
303      CoinBigIndex startRow=rowStart[iRow];
304      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
305      work[iRow] = diagonalSlack[iRow]+delta2;
306      for (CoinBigIndex k=startRow;k<endRow;k++) {
307        int iColumn=column[k];
308        if (!whichDense_||!whichDense_[iColumn]) {
309          CoinBigIndex start=columnStart[iColumn];
310          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
311          double multiplier = diagonal[iColumn]*elementByRow[k];
312          for (CoinBigIndex j=start;j<end;j++) {
313            int jRow=row[j];
314            if (jRow>=iRow&&!rowsDropped_[jRow]) {
315              double value=element[j]*multiplier;
316              work[jRow] += value;
317            }
318          }
319        }
320      }
321      int j;
322      for (j=0;j<number;j++) {
323        int jRow = which[j];
324        put[j]=work[jRow];
325        work[jRow]=0.0;
326      }
327    } else {
328      // dropped
329      int j;
330      for (j=1;j<number;j++) {
331        put[j]=0.0;
332      }
333      put[0]=1.0;
334    }
335  }
336  //check sizes
337  double largest2=maximumAbsElement(sparseFactor_,sizeFactor_);
338  largest2*=1.0e-20;
339  largest = CoinMin(largest2,1.0e-11);
340  int numberDroppedBefore=0;
341  for (iRow=0;iRow<numberRows_;iRow++) {
342    int dropped=rowsDropped_[iRow];
343    // Move to int array
344    rowsDropped[iRow]=dropped;
345    if (!dropped) {
346      CoinBigIndex start = choleskyStart_[iRow];
347      double diagonal = sparseFactor_[start];
348      if (diagonal>largest2) {
349        sparseFactor_[start]=CoinMax(diagonal,1.0e-10);
350      } else {
351        sparseFactor_[start]=CoinMax(diagonal,1.0e-10);
352        rowsDropped[iRow]=2;
353        numberDroppedBefore++;
354      } 
355    } 
356  } 
357  delete [] work;
358  cholmod_sparse A ;
359  A.nrow=numberRows_;
360  A.ncol=numberRows_;
361  A.nzmax=choleskyStart_[numberRows_];
362  A.p = choleskyStart_;
363  A.i = choleskyRow_;
364  A.x=sparseFactor_;
365  A.stype=-1;
366  A.itype=CHOLMOD_INT;
367  A.xtype=CHOLMOD_REAL;
368  A.dtype=CHOLMOD_DOUBLE;
369  A.sorted=1;
370  A.packed=1;
371  cholmod_factorize (&A, L_, &c_) ;                 /* factorize */
372  choleskyCondition_=1.0;
373  bool cleanCholesky;
374  if (model_->numberIterations()<2000) 
375    cleanCholesky=true;
376  else 
377    cleanCholesky=false;
378  if (cleanCholesky) {
379    //drop fresh makes some formADAT easier
380    //int oldDropped=numberRowsDropped_;
381    if (newDropped||numberRowsDropped_) {
382      //std::cout <<"Rank "<<numberRows_-newDropped<<" ( "<<
383      //  newDropped<<" dropped)";
384      //if (newDropped>oldDropped)
385      //std::cout<<" ( "<<newDropped-oldDropped<<" dropped this time)";
386      //std::cout<<std::endl;
387      newDropped=0;
388      for (int i=0;i<numberRows_;i++) {
389        char dropped = rowsDropped[i];
390        rowsDropped_[i]=dropped;
391        if (dropped==2) {
392          //dropped this time
393          rowsDropped[newDropped++]=i;
394          rowsDropped_[i]=0;
395        } 
396      } 
397      numberRowsDropped_=newDropped;
398      newDropped=-(2+newDropped);
399    } 
400  } else {
401    if (newDropped) {
402      newDropped=0;
403      for (int i=0;i<numberRows_;i++) {
404        char dropped = rowsDropped[i];
405        rowsDropped_[i]=dropped;
406        if (dropped==2) {
407          //dropped this time
408          rowsDropped[newDropped++]=i;
409          rowsDropped_[i]=1;
410        } 
411      } 
412    } 
413    numberRowsDropped_+=newDropped;
414    if (numberRowsDropped_&&0) {
415      std::cout <<"Rank "<<numberRows_-numberRowsDropped_<<" ( "<<
416          numberRowsDropped_<<" dropped)";
417      if (newDropped) {
418        std::cout<<" ( "<<newDropped<<" dropped this time)";
419      } 
420      std::cout<<std::endl;
421    } 
422  }
423  status_=0;
424  return newDropped;
425}
426/* Uses factorization to solve. */
427void 
428ClpCholeskyUfl::solve (double * region) 
429{
430  cholmod_dense *x, *b;
431  b = cholmod_allocate_dense (numberRows_, 1, numberRows_, CHOLMOD_REAL, &c_) ; 
432  CoinMemcpyN(region,numberRows_,(double *) b->x);
433  x = cholmod_solve (CHOLMOD_A, L_, b, &c_) ;
434  CoinMemcpyN((double *) x->x,numberRows_,region);
435  cholmod_free_dense (&x, &c_) ;
436  cholmod_free_dense (&b, &c_) ;
437}
438#endif
439#endif
Note: See TracBrowser for help on using the repository browser.