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

Last change on this file since 1355 was 423, checked in by forrest, 16 years ago

various Cholesky

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.9 KB
Line 
1#ifdef TAUCS_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 "ClpCholeskyTaucs.hpp"
13#include "ClpMessage.hpp"
14
15//#############################################################################
16// Constructors / Destructor / Assignment
17//#############################################################################
18
19//-------------------------------------------------------------------
20// Default Constructor
21//-------------------------------------------------------------------
22ClpCholeskyTaucs::ClpCholeskyTaucs () 
23  : ClpCholeskyBase(),
24    matrix_(NULL),
25    factorization_(NULL),
26    sparseFactorT_(NULL),
27    choleskyStartT_(NULL),
28    choleskyRowT_(NULL),
29    sizeFactorT_(0),
30    rowCopyT_(NULL)
31{
32  type_=13;
33}
34
35//-------------------------------------------------------------------
36// Copy constructor
37//-------------------------------------------------------------------
38ClpCholeskyTaucs::ClpCholeskyTaucs (const ClpCholeskyTaucs & rhs) 
39: ClpCholeskyBase(rhs)
40{
41  type_=rhs.type_;
42  // For Taucs stuff is done by malloc
43  matrix_ = rhs.matrix_;
44  sizeFactorT_=rhs.sizeFactorT_;
45  if (matrix_) {
46    choleskyStartT_ = (int *) malloc((numberRows_+1)*sizeof(int));
47    memcpy(choleskyStartT_,rhs.choleskyStartT_,(numberRows_+1)*sizeof(int));
48    choleskyRowT_ = (int *) malloc(sizeFactorT_*sizeof(int));
49    memcpy(choleskyRowT_,rhs.choleskyRowT_,sizeFactorT_*sizeof(int));
50    sparseFactorT_ = (double *) malloc(sizeFactorT_*sizeof(double));
51    memcpy(sparseFactorT_,rhs.sparseFactorT_,sizeFactorT_*sizeof(double));
52    matrix_->colptr = choleskyStartT_;
53    matrix_->rowind = choleskyRowT_;
54    matrix_->values.d = sparseFactorT_;
55  } else {
56    sparseFactorT_=NULL;
57    choleskyStartT_=NULL;
58    choleskyRowT_=NULL;
59  }
60  factorization_=NULL,
61  rowCopyT_ = rhs.rowCopyT_->clone();
62}
63
64
65//-------------------------------------------------------------------
66// Destructor
67//-------------------------------------------------------------------
68ClpCholeskyTaucs::~ClpCholeskyTaucs ()
69{
70  taucs_ccs_free(matrix_);
71  if (factorization_)
72    taucs_supernodal_factor_free(factorization_);
73  delete rowCopyT_;
74}
75
76//----------------------------------------------------------------
77// Assignment operator
78//-------------------------------------------------------------------
79ClpCholeskyTaucs &
80ClpCholeskyTaucs::operator=(const ClpCholeskyTaucs& rhs)
81{
82  if (this != &rhs) {
83    ClpCholeskyBase::operator=(rhs);
84    taucs_ccs_free(matrix_);
85    if (factorization_)
86      taucs_supernodal_factor_free(factorization_);
87    factorization_=NULL;
88    sizeFactorT_=rhs.sizeFactorT_;
89    matrix_ = rhs.matrix_;
90    if (matrix_) {
91      choleskyStartT_ = (int *) malloc((numberRows_+1)*sizeof(int));
92      memcpy(choleskyStartT_,rhs.choleskyStartT_,(numberRows_+1)*sizeof(int));
93      choleskyRowT_ = (int *) malloc(sizeFactorT_*sizeof(int));
94      memcpy(choleskyRowT_,rhs.choleskyRowT_,sizeFactorT_*sizeof(int));
95      sparseFactorT_ = (double *) malloc(sizeFactorT_*sizeof(double));
96      memcpy(sparseFactorT_,rhs.sparseFactorT_,sizeFactorT_*sizeof(double));
97      matrix_->colptr = choleskyStartT_;
98      matrix_->rowind = choleskyRowT_;
99      matrix_->values.d = sparseFactorT_;
100    } else {
101      sparseFactorT_=NULL;
102      choleskyStartT_=NULL;
103      choleskyRowT_=NULL;
104    }
105    delete rowCopyT_;
106    rowCopyT_ = rhs.rowCopyT_->clone();
107  }
108  return *this;
109}
110//-------------------------------------------------------------------
111// Clone
112//-------------------------------------------------------------------
113ClpCholeskyBase * ClpCholeskyTaucs::clone() const
114{
115  return new ClpCholeskyTaucs(*this);
116}
117/* Orders rows and saves pointer to matrix.and model */
118int 
119ClpCholeskyTaucs::order(ClpInterior * model) 
120{
121  numberRows_ = model->numberRows();
122  rowsDropped_ = new char [numberRows_];
123  memset(rowsDropped_,0,numberRows_);
124  numberRowsDropped_=0;
125  model_=model;
126  rowCopyT_ = model->clpMatrix()->reverseOrderedCopy();
127  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
128  const int * columnLength = model_->clpMatrix()->getVectorLengths();
129  const int * row = model_->clpMatrix()->getIndices();
130  const CoinBigIndex * rowStart = rowCopyT_->getVectorStarts();
131  const int * rowLength = rowCopyT_->getVectorLengths();
132  const int * column = rowCopyT_->getIndices();
133  // We need two arrays for counts
134  int * which = new int [numberRows_];
135  int * used = new int[numberRows_];
136  CoinZeroN(used,numberRows_);
137  int iRow;
138  sizeFactorT_=0;
139  for (iRow=0;iRow<numberRows_;iRow++) {
140    int number=1;
141    // make sure diagonal exists
142    which[0]=iRow;
143    used[iRow]=1;
144    if (!rowsDropped_[iRow]) {
145      CoinBigIndex startRow=rowStart[iRow];
146      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
147      for (CoinBigIndex k=startRow;k<endRow;k++) {
148        int iColumn=column[k];
149        CoinBigIndex start=columnStart[iColumn];
150        CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
151        for (CoinBigIndex j=start;j<end;j++) {
152          int jRow=row[j];
153          if (jRow>=iRow&&!rowsDropped_[jRow]) {
154            if (!used[jRow]) {
155              used[jRow]=1;
156              which[number++]=jRow;
157            }
158          }
159        }
160      }
161      sizeFactorT_ += number;
162      int j;
163      for (j=0;j<number;j++)
164        used[which[j]]=0;
165    }
166  }
167  delete [] which;
168  // Now we have size - create arrays and fill in
169  matrix_ = taucs_ccs_create(numberRows_,numberRows_,sizeFactorT_,
170                             TAUCS_DOUBLE|TAUCS_SYMMETRIC|TAUCS_LOWER); 
171  if (!matrix_) 
172    return 1;
173  // Space for starts
174  choleskyStartT_ = matrix_->colptr;
175  choleskyRowT_ = matrix_->rowind;
176  sparseFactorT_ = matrix_->values.d;
177  sizeFactorT_=0;
178  which = choleskyRowT_;
179  for (iRow=0;iRow<numberRows_;iRow++) {
180    int number=1;
181    // make sure diagonal exists
182    which[0]=iRow;
183    used[iRow]=1;
184    choleskyStartT_[iRow]=sizeFactorT_;
185    if (!rowsDropped_[iRow]) {
186      CoinBigIndex startRow=rowStart[iRow];
187      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
188      for (CoinBigIndex k=startRow;k<endRow;k++) {
189        int iColumn=column[k];
190        CoinBigIndex start=columnStart[iColumn];
191        CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
192        for (CoinBigIndex j=start;j<end;j++) {
193          int jRow=row[j];
194          if (jRow>=iRow&&!rowsDropped_[jRow]) {
195            if (!used[jRow]) {
196              used[jRow]=1;
197              which[number++]=jRow;
198            }
199          }
200        }
201      }
202      sizeFactorT_ += number;
203      int j;
204      for (j=0;j<number;j++)
205        used[which[j]]=0;
206      // Sort
207      std::sort(which,which+number);
208      // move which on
209      which += number;
210    }
211  }
212  choleskyStartT_[numberRows_]=sizeFactorT_;
213  delete [] used;
214  permuteInverse_ = new int [numberRows_];
215  permute_ = new int[numberRows_];
216  int * perm, *invp;
217  // There seem to be bugs in ordering if model too small
218  if (numberRows_>10)
219    taucs_ccs_order(matrix_,&perm,&invp,(const char *) "genmmd");
220  else
221    taucs_ccs_order(matrix_,&perm,&invp,(const char *) "identity");
222  memcpy(permuteInverse_,perm,numberRows_*sizeof(int));
223  free(perm);
224  memcpy(permute_,invp,numberRows_*sizeof(int));
225  free(invp);
226  // need to permute
227  taucs_ccs_matrix * permuted = taucs_ccs_permute_symmetrically(matrix_,permuteInverse_,permute_);
228  // symbolic
229  factorization_ = taucs_ccs_factor_llt_symbolic(permuted);
230  taucs_ccs_free(permuted);
231  return factorization_ ? 0 :1;
232}
233/* Does Symbolic factorization given permutation.
234   This is called immediately after order.  If user provides this then
235   user must provide factorize and solve.  Otherwise the default factorization is used
236   returns non-zero if not enough memory */
237int 
238ClpCholeskyTaucs::symbolic()
239{
240  return 0;
241}
242/* Factorize - filling in rowsDropped and returning number dropped */
243int 
244ClpCholeskyTaucs::factorize(const double * diagonal, int * rowsDropped) 
245{
246  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
247  const int * columnLength = model_->clpMatrix()->getVectorLengths();
248  const int * row = model_->clpMatrix()->getIndices();
249  const double * element = model_->clpMatrix()->getElements();
250  const CoinBigIndex * rowStart = rowCopyT_->getVectorStarts();
251  const int * rowLength = rowCopyT_->getVectorLengths();
252  const int * column = rowCopyT_->getIndices();
253  const double * elementByRow = rowCopyT_->getElements();
254  int numberColumns=model_->clpMatrix()->getNumCols();
255  int iRow;
256  double * work = new double[numberRows_];
257  CoinZeroN(work,numberRows_);
258  const double * diagonalSlack = diagonal + numberColumns;
259  int newDropped=0;
260  double largest;
261  //perturbation
262  double perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
263  perturbation=perturbation*perturbation;
264  if (perturbation>1.0) {
265    //if (model_->model()->logLevel()&4)
266      std::cout <<"large perturbation "<<perturbation<<std::endl;
267    perturbation=sqrt(perturbation);;
268    perturbation=1.0;
269  } 
270  for (iRow=0;iRow<numberRows_;iRow++) {
271    double * put = sparseFactorT_+choleskyStartT_[iRow];
272    int * which = choleskyRowT_+choleskyStartT_[iRow];
273    int number = choleskyStartT_[iRow+1]-choleskyStartT_[iRow];
274    if (!rowLength[iRow])
275      rowsDropped_[iRow]=1;
276    if (!rowsDropped_[iRow]) {
277      CoinBigIndex startRow=rowStart[iRow];
278      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
279      work[iRow] = diagonalSlack[iRow];
280      for (CoinBigIndex k=startRow;k<endRow;k++) {
281        int iColumn=column[k];
282        CoinBigIndex start=columnStart[iColumn];
283        CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
284        double multiplier = diagonal[iColumn]*elementByRow[k];
285        for (CoinBigIndex j=start;j<end;j++) {
286          int jRow=row[j];
287          if (jRow>=iRow&&!rowsDropped_[jRow]) {
288            double value=element[j]*multiplier;
289            work[jRow] += value;
290          }
291        }
292      }
293      int j;
294      for (j=0;j<number;j++) {
295        int jRow = which[j];
296        put[j]=work[jRow];
297        work[jRow]=0.0;
298      }
299    } else {
300      // dropped
301      int j;
302      for (j=1;j<number;j++) {
303        put[j]=0.0;
304      }
305      put[0]=1.0;
306    }
307  }
308  //check sizes
309  double largest2=maximumAbsElement(sparseFactorT_,sizeFactorT_);
310  largest2*=1.0e-19;
311  largest = CoinMin(largest2,1.0e-11);
312  int numberDroppedBefore=0;
313  for (iRow=0;iRow<numberRows_;iRow++) {
314    int dropped=rowsDropped_[iRow];
315    // Move to int array
316    rowsDropped[iRow]=dropped;
317    if (!dropped) {
318      CoinBigIndex start = choleskyStartT_[iRow];
319      double diagonal = sparseFactorT_[start];
320      if (diagonal>largest2) {
321        sparseFactorT_[start]=diagonal+perturbation;
322      } else {
323        sparseFactorT_[start]=diagonal+perturbation;
324        rowsDropped[iRow]=2;
325        numberDroppedBefore++;
326      } 
327    } 
328  }
329  taucs_supernodal_factor_free_numeric(factorization_);
330  // need to permute
331  taucs_ccs_matrix * permuted = taucs_ccs_permute_symmetrically(matrix_,permuteInverse_,permute_);
332  int rCode=taucs_ccs_factor_llt_numeric(permuted,factorization_);
333  taucs_ccs_free(permuted);
334  if (rCode)
335    printf("return code of %d from factor\n",rCode);
336  delete [] work;
337  choleskyCondition_=1.0;
338  bool cleanCholesky;
339  if (model_->numberIterations()<200) 
340    cleanCholesky=true;
341  else 
342    cleanCholesky=false;
343  /*
344    How do I find out where 1.0e100's are in cholesky?
345  */
346  if (cleanCholesky) {
347    //drop fresh makes some formADAT easier
348    int oldDropped=numberRowsDropped_;
349    if (newDropped||numberRowsDropped_) {
350      std::cout <<"Rank "<<numberRows_-newDropped<<" ( "<<
351          newDropped<<" dropped)";
352      if (newDropped>oldDropped) 
353        std::cout<<" ( "<<newDropped-oldDropped<<" dropped this time)";
354      std::cout<<std::endl;
355      newDropped=0;
356      for (int i=0;i<numberRows_;i++) {
357        char dropped = rowsDropped[i];
358        rowsDropped_[i]=dropped;
359        if (dropped==2) {
360          //dropped this time
361          rowsDropped[newDropped++]=i;
362          rowsDropped_[i]=0;
363        } 
364      } 
365      numberRowsDropped_=newDropped;
366      newDropped=-(1+newDropped);
367    } 
368  } else {
369    if (newDropped) {
370      newDropped=0;
371      for (int i=0;i<numberRows_;i++) {
372        char dropped = rowsDropped[i];
373        int oldDropped = rowsDropped_[i];;
374        rowsDropped_[i]=dropped;
375        if (dropped==2) {
376          assert (!oldDropped);
377          //dropped this time
378          rowsDropped[newDropped++]=i;
379          rowsDropped_[i]=1;
380        } 
381      } 
382    } 
383    numberRowsDropped_+=newDropped;
384    if (numberRowsDropped_) {
385      std::cout <<"Rank "<<numberRows_-numberRowsDropped_<<" ( "<<
386          numberRowsDropped_<<" dropped)";
387      if (newDropped) {
388        std::cout<<" ( "<<newDropped<<" dropped this time)";
389      } 
390      std::cout<<std::endl;
391    } 
392  } 
393  status_=0;
394  return newDropped;
395}
396/* Uses factorization to solve. */
397void 
398ClpCholeskyTaucs::solve (double * region) 
399{
400  double * in = new double[numberRows_];
401  double * out = new double[numberRows_];
402  taucs_vec_permute(numberRows_,TAUCS_DOUBLE,region,in,permuteInverse_);
403  int rCode=taucs_supernodal_solve_llt(factorization_,out,in);
404  if (rCode)
405    printf("return code of %d from solve\n",rCode);
406  taucs_vec_permute(numberRows_,TAUCS_DOUBLE,out,region,permute_);
407  delete [] out;
408  delete [] in;
409}
410#endif
Note: See TracBrowser for help on using the repository browser.