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

Last change on this file since 800 was 800, checked in by ladanyi, 14 years ago

finishing conversion to svn

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