source: trunk/ClpCholeskyUfl.cpp @ 719

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

minor stuff

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