source: stable/1.15/Clp/src/ClpCholeskyUfl.cpp @ 1989

Last change on this file since 1989 was 1723, checked in by forrest, 9 years ago

out some printf statements

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