00001
00002
00003
00004
00005
00006 #include "CoinPragma.hpp"
00007 #include "CoinHelperFunctions.hpp"
00008
00009 #include "ClpInterior.hpp"
00010 #include "ClpCholeskyDense.hpp"
00011 #include "ClpMessage.hpp"
00012
00013
00014
00015
00016
00017
00018
00019
00020 ClpCholeskyDense::ClpCholeskyDense ()
00021 : ClpCholeskyBase(),
00022 work_(NULL),
00023 rowCopy_(NULL)
00024 {
00025 type_=11;;
00026 }
00027
00028
00029
00030
00031 ClpCholeskyDense::ClpCholeskyDense (const ClpCholeskyDense & rhs)
00032 : ClpCholeskyBase(rhs)
00033 {
00034 type_=rhs.type_;
00035 work_ = ClpCopyOfArray(rhs.work_,numberRows_*numberRows_);
00036 rowCopy_ = rhs.rowCopy_->clone();
00037 }
00038
00039
00040
00041
00042
00043 ClpCholeskyDense::~ClpCholeskyDense ()
00044 {
00045 delete [] work_;
00046 delete rowCopy_;
00047 }
00048
00049
00050
00051
00052 ClpCholeskyDense &
00053 ClpCholeskyDense::operator=(const ClpCholeskyDense& rhs)
00054 {
00055 if (this != &rhs) {
00056 ClpCholeskyBase::operator=(rhs);
00057 delete [] work_;
00058 work_ = ClpCopyOfArray(rhs.work_,numberRows_*numberRows_);
00059 delete rowCopy_;
00060 rowCopy_ = rhs.rowCopy_->clone();
00061 }
00062 return *this;
00063 }
00064
00065
00066
00067 ClpCholeskyBase * ClpCholeskyDense::clone() const
00068 {
00069 return new ClpCholeskyDense(*this);
00070 }
00071
00072 int
00073 ClpCholeskyDense::order(ClpInterior * model)
00074 {
00075 numberRows_ = model->numberRows();
00076 work_ = new double [numberRows_*numberRows_];
00077 rowsDropped_ = new char [numberRows_];
00078 memset(rowsDropped_,0,numberRows_);
00079 numberRowsDropped_=0;
00080 model_=model;
00081 rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
00082 return 0;
00083 }
00084
00085
00086 int
00087 ClpCholeskyDense::factorize(const double * diagonal, int * rowsDropped)
00088 {
00089 int iColumn;
00090 const CoinBigIndex * columnStart = model_->matrix()->getVectorStarts();
00091 const int * columnLength = model_->matrix()->getVectorLengths();
00092 const int * row = model_->matrix()->getIndices();
00093 const double * element = model_->matrix()->getElements();
00094 const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
00095 const int * rowLength = rowCopy_->getVectorLengths();
00096 const int * column = rowCopy_->getIndices();
00097 const double * elementByRow = rowCopy_->getElements();
00098 int numberColumns=model_->matrix()->getNumCols();
00099 CoinZeroN(work_,numberRows_*numberRows_);
00100 int iRow;
00101 double * work = work_;
00102 const double * diagonalSlack = diagonal + numberColumns;
00103 for (iRow=0;iRow<numberRows_;iRow++) {
00104 if (!rowsDropped_[iRow]) {
00105 CoinBigIndex startRow=rowStart[iRow];
00106 CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
00107 work[iRow] = diagonalSlack[iRow];
00108 for (CoinBigIndex k=startRow;k<endRow;k++) {
00109 int iColumn=column[k];
00110 CoinBigIndex start=columnStart[iColumn];
00111 CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
00112 double multiplier = diagonal[iColumn]*elementByRow[k];
00113 for (CoinBigIndex j=start;j<end;j++) {
00114 int jRow=row[j];
00115 if (jRow>=iRow&&!rowsDropped_[jRow]) {
00116 double value=element[j]*multiplier;
00117 work[jRow] += value;
00118 }
00119 }
00120 }
00121 }
00122 work += numberRows_;
00123 }
00124 work = work_;
00125 #ifdef CLP_DEBUG
00126 double * save = NULL;
00127 if (numberRows_<20) {
00128 save = new double [ numberRows_*numberRows_];
00129 memcpy(save,work,numberRows_*numberRows_*sizeof(double));
00130 }
00131 #endif
00132 int numberDropped=0;
00133 for (iColumn=0;iColumn<numberRows_;iColumn++) {
00134 int iRow;
00135 if (!rowsDropped_[iColumn]) {
00136 double diagonalValue = work[iColumn];
00137 if (diagonalValue>pivotTolerance_) {
00138 diagonalValue = sqrt(diagonalValue);
00139 work[iColumn]=diagonalValue;
00140 for (iRow=iColumn+1;iRow<numberRows_;iRow++)
00141 work[iRow] /= diagonalValue;
00142 double * work2 = work;
00143 for (int jColumn=iColumn+1;jColumn<numberRows_;jColumn++) {
00144 work2 += numberRows_;
00145 double value = work[jColumn];
00146 for (iRow=jColumn;iRow<numberRows_;iRow++)
00147 work2[iRow] -= value*work[iRow];
00148 }
00149 } else {
00150
00151 rowsDropped_[iColumn]=-1;
00152 rowsDropped[numberDropped++]=iColumn;
00153 numberRowsDropped_++;
00154
00155 for (iRow=0;iRow<iColumn;iRow++)
00156 work_[iColumn+iRow*numberRows_]=0.0;
00157 for (iRow=iColumn;iRow<numberRows_;iRow++)
00158 work[iRow]=0.0;
00159 }
00160 } else {
00161
00162 for (iRow=0;iRow<iColumn;iRow++)
00163 work_[iColumn+iRow*numberRows_]=0.0;
00164 for (iRow=iColumn;iRow<numberRows_;iRow++)
00165 work[iRow]=0.0;
00166 }
00167 work += numberRows_;
00168 }
00169 #ifdef CLP_DEBUG
00170 if (save) {
00171 double * array = new double [numberRows_];
00172 int iColumn;
00173 for (iColumn=0;iColumn<numberRows_;iColumn++) {
00174 double * s = save;
00175 int iRow;
00176 for (iRow=0;iRow<iColumn;iRow++) {
00177 array[iRow]=s[iColumn];
00178 s += numberRows_;
00179 }
00180 for (iRow=iColumn;iRow<numberRows_;iRow++) {
00181 array[iRow]=s[iRow];
00182 }
00183 solve(array);
00184 for (iRow=0;iRow<numberRows_;iRow++) {
00185 if (iRow!=iColumn)
00186 assert (fabs(array[iRow])<1.0e-7);
00187 else
00188 assert (fabs(array[iRow]-1.0)<1.0e-7);
00189 }
00190 }
00191 delete [] array;
00192 delete [] save;
00193 }
00194 #endif
00195 return numberDropped;
00196 }
00197
00198 void
00199 ClpCholeskyDense::solve (double * region) const
00200 {
00201 int iColumn;
00202 for (iColumn=0;iColumn<numberRows_;iColumn++) {
00203 if (!rowsDropped_[iColumn]) {
00204 double value = region[iColumn];
00205 for (int iRow=0;iRow<iColumn;iRow++)
00206 value -= region[iRow]*work_[iColumn+iRow*numberRows_];
00207 for (int iRow=0;iRow<iColumn;iRow++)
00208 if (rowsDropped_[iRow])
00209 assert(!work_[iColumn+iRow*numberRows_]||!region[iRow]);
00210 region[iColumn]=value/work_[iColumn+iColumn*numberRows_];
00211 } else {
00212 region[iColumn]=0.0;
00213 }
00214 }
00215 double * work = work_ + numberRows_*numberRows_;
00216 for (iColumn=numberRows_-1;iColumn>=0;iColumn--) {
00217 work -= numberRows_;
00218 if (!rowsDropped_[iColumn]) {
00219 double value = region[iColumn];
00220 for (int iRow=iColumn+1;iRow<numberRows_;iRow++)
00221 value -= region[iRow]*work[iRow];
00222 for (int iRow=iColumn+1;iRow<numberRows_;iRow++)
00223 if (rowsDropped_[iRow])
00224 assert(!work[iRow]||!region[iRow]);
00225 region[iColumn]=value/work[iColumn];
00226 } else {
00227 region[iColumn]=0.0;
00228 }
00229 }
00230 }