source: trunk/ClpCholeskyDense.cpp @ 284

Last change on this file since 284 was 284, checked in by jpfasano, 16 years ago

Modified to compile with MS Visual Studio Version 6

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.9 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5
6#include "CoinPragma.hpp"
7#include "CoinHelperFunctions.hpp"
8
9#include "ClpInterior.hpp"
10#include "ClpCholeskyDense.hpp"
11#include "ClpMessage.hpp"
12
13//#############################################################################
14// Constructors / Destructor / Assignment
15//#############################################################################
16
17//-------------------------------------------------------------------
18// Default Constructor
19//-------------------------------------------------------------------
20ClpCholeskyDense::ClpCholeskyDense () 
21  : ClpCholeskyBase(),
22    work_(NULL),
23    rowCopy_(NULL)
24{
25  type_=11;;
26}
27
28//-------------------------------------------------------------------
29// Copy constructor
30//-------------------------------------------------------------------
31ClpCholeskyDense::ClpCholeskyDense (const ClpCholeskyDense & rhs) 
32: ClpCholeskyBase(rhs)
33{
34  type_=rhs.type_;
35  work_ = ClpCopyOfArray(rhs.work_,numberRows_*numberRows_);
36  rowCopy_ = rhs.rowCopy_->clone();
37}
38
39
40//-------------------------------------------------------------------
41// Destructor
42//-------------------------------------------------------------------
43ClpCholeskyDense::~ClpCholeskyDense ()
44{
45  delete [] work_;
46  delete rowCopy_;
47}
48
49//----------------------------------------------------------------
50// Assignment operator
51//-------------------------------------------------------------------
52ClpCholeskyDense &
53ClpCholeskyDense::operator=(const ClpCholeskyDense& rhs)
54{
55  if (this != &rhs) {
56    ClpCholeskyBase::operator=(rhs);
57    delete [] work_;
58    work_ = ClpCopyOfArray(rhs.work_,numberRows_*numberRows_);
59    delete rowCopy_;
60    rowCopy_ = rhs.rowCopy_->clone();
61  }
62  return *this;
63}
64//-------------------------------------------------------------------
65// Clone
66//-------------------------------------------------------------------
67ClpCholeskyBase * ClpCholeskyDense::clone() const
68{
69  return new ClpCholeskyDense(*this);
70}
71/* Orders rows and saves pointer to matrix.and model */
72int 
73ClpCholeskyDense::order(ClpInterior * model) 
74{
75  numberRows_ = model->numberRows();
76  work_ = new double [numberRows_*numberRows_];
77  rowsDropped_ = new char [numberRows_];
78  memset(rowsDropped_,0,numberRows_);
79  numberRowsDropped_=0;
80  model_=model;
81  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
82  return 0;
83}
84//#define CLP_DEBUG
85/* Factorize - filling in rowsDropped and returning number dropped */
86int 
87ClpCholeskyDense::factorize(const double * diagonal, int * rowsDropped) 
88{
89  int iColumn;
90  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
91  const int * columnLength = model_->clpMatrix()->getVectorLengths();
92  const int * row = model_->clpMatrix()->getIndices();
93  const double * element = model_->clpMatrix()->getElements();
94  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
95  const int * rowLength = rowCopy_->getVectorLengths();
96  const int * column = rowCopy_->getIndices();
97  const double * elementByRow = rowCopy_->getElements();
98  int numberColumns=model_->clpMatrix()->getNumCols();
99  CoinZeroN(work_,numberRows_*numberRows_);
100  int iRow;
101  double * work = work_;
102  const double * diagonalSlack = diagonal + numberColumns;
103  for (iRow=0;iRow<numberRows_;iRow++) {
104    if (!rowsDropped_[iRow]) {
105      CoinBigIndex startRow=rowStart[iRow];
106      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
107      work[iRow] = diagonalSlack[iRow];
108      for (CoinBigIndex k=startRow;k<endRow;k++) {
109        int iColumn=column[k];
110        CoinBigIndex start=columnStart[iColumn];
111        CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
112        double multiplier = diagonal[iColumn]*elementByRow[k];
113        for (CoinBigIndex j=start;j<end;j++) {
114          int jRow=row[j];
115          if (jRow>=iRow&&!rowsDropped_[jRow]) {
116            double value=element[j]*multiplier;
117            work[jRow] += value;
118          }
119        }
120      } 
121    }
122    work += numberRows_;
123  }
124  work = work_;
125  // Get rid of any odd copies
126  // model_->clpMatrix()->releasePackedMatrix();
127#ifdef CLP_DEBUG
128  double * save = NULL;
129  if (numberRows_<20) {
130    save = new double [ numberRows_*numberRows_];
131    memcpy(save,work,numberRows_*numberRows_*sizeof(double));
132  }
133#endif
134  int numberDropped=0;
135  for (iColumn=0;iColumn<numberRows_;iColumn++) {
136    int iRow;
137    if (!rowsDropped_[iColumn]) {
138      double diagonalValue = work[iColumn];
139      if (diagonalValue>pivotTolerance_) {
140        diagonalValue = sqrt(diagonalValue);
141        work[iColumn]=diagonalValue;
142        for (iRow=iColumn+1;iRow<numberRows_;iRow++)
143          work[iRow] /= diagonalValue;
144        double * work2 = work;
145        for (int jColumn=iColumn+1;jColumn<numberRows_;jColumn++) {
146          work2 += numberRows_;
147          double value = work[jColumn];
148          for (iRow=jColumn;iRow<numberRows_;iRow++)
149          work2[iRow] -= value*work[iRow];
150        }
151      } else {
152        // drop column
153        rowsDropped_[iColumn]=-1;
154        rowsDropped[numberDropped++]=iColumn;
155        numberRowsDropped_++;
156        // clean up as this is a debug version
157        for (iRow=0;iRow<iColumn;iRow++)
158          work_[iColumn+iRow*numberRows_]=0.0;
159        for (iRow=iColumn;iRow<numberRows_;iRow++)
160          work[iRow]=0.0;
161      }
162    } else {
163      // clean up as this is a debug version
164      for (iRow=0;iRow<iColumn;iRow++)
165        work_[iColumn+iRow*numberRows_]=0.0;
166      for (iRow=iColumn;iRow<numberRows_;iRow++)
167        work[iRow]=0.0;
168    }
169    work += numberRows_; // move on pointer
170  }
171#ifdef CLP_DEBUG
172  if (save) {
173    double * array = new double [numberRows_];
174    int iColumn;
175    for (iColumn=0;iColumn<numberRows_;iColumn++) {
176      double * s = save;
177      int iRow;
178      for (iRow=0;iRow<iColumn;iRow++) {
179        array[iRow]=s[iColumn];
180        s += numberRows_;
181      }
182      for (iRow=iColumn;iRow<numberRows_;iRow++) {
183        array[iRow]=s[iRow];
184      }
185      solve(array);
186      for (iRow=0;iRow<numberRows_;iRow++) {
187        if (iRow!=iColumn)
188          assert (fabs(array[iRow])<1.0e-7);
189        else
190          assert (fabs(array[iRow]-1.0)<1.0e-7);
191      }
192    }
193    delete [] array;
194    delete [] save;
195  }
196#endif
197  return numberDropped;
198}
199/* Uses factorization to solve. */
200void 
201ClpCholeskyDense::solve (double * region) 
202{
203  int iColumn;
204  for (iColumn=0;iColumn<numberRows_;iColumn++) {
205    if (!rowsDropped_[iColumn]) {
206      double value = region[iColumn];
207          int iRow;
208      for (iRow=0;iRow<iColumn;iRow++)
209        value -= region[iRow]*work_[iColumn+iRow*numberRows_];
210      for (iRow=0;iRow<iColumn;iRow++)
211        if (rowsDropped_[iRow])
212          assert(!work_[iColumn+iRow*numberRows_]||!region[iRow]);
213      region[iColumn]=value/work_[iColumn+iColumn*numberRows_];
214    } else {
215      region[iColumn]=0.0;
216    }
217  }
218  double * work = work_ + numberRows_*numberRows_;
219  for (iColumn=numberRows_-1;iColumn>=0;iColumn--) {
220    work -= numberRows_;
221    if (!rowsDropped_[iColumn]) {
222      double value = region[iColumn];
223      for (int iRow=iColumn+1;iRow<numberRows_;iRow++)
224        value -= region[iRow]*work[iRow];
225      for (iRow=iColumn+1;iRow<numberRows_;iRow++)
226        if (rowsDropped_[iRow])
227          assert(!work[iRow]||!region[iRow]);
228      region[iColumn]=value/work[iColumn];
229    } else {
230      region[iColumn]=0.0;
231    }
232  }
233}
Note: See TracBrowser for help on using the repository browser.