source: trunk/Clp/src/ClpCholeskyMumps.cpp @ 1665

Last change on this file since 1665 was 1665, checked in by lou, 9 years ago

Add EPL license notice in src.

  • Property svn:keywords set to Id
File size: 16.0 KB
Line 
1/* $Id: ClpCholeskyMumps.cpp 1665 2011-01-04 17:55:54Z lou $ */
2#ifdef MUMPS_BARRIER
3// Copyright (C) 2004, International Business Machines
4// Corporation and others.  All Rights Reserved.
5// This code is licensed under the terms of the Eclipse Public License (EPL).
6
7
8#include "CoinPragma.hpp"
9#include "ClpCholeskyMumps.hpp"
10extern "C"
11{
12     //#include "mpi.h"
13#include "/home/jjforre/cbc-trunk/ThirdParty/Mumps/MUMPS/libseq/mpi.h"
14}
15#define JOB_INIT -1
16#define JOB_END -2
17#define USE_COMM_WORLD -987654
18#include "ClpMessage.hpp"
19#include "ClpInterior.hpp"
20#include "CoinHelperFunctions.hpp"
21#include "ClpHelperFunctions.hpp"
22//#############################################################################
23// Constructors / Destructor / Assignment
24//#############################################################################
25
26//-------------------------------------------------------------------
27// Default Constructor
28//-------------------------------------------------------------------
29ClpCholeskyMumps::ClpCholeskyMumps (int denseThreshold)
30     : ClpCholeskyBase(denseThreshold)
31{
32     type_ = 16;
33     mumps_.n = 0;
34     mumps_.nz = 0;
35     mumps_.a = NULL;
36     mumps_.jcn = NULL;
37     mumps_.irn = NULL;
38     mumps_.job = JOB_INIT;//initialize mumps
39     mumps_.par = 1;//working host for sequential version
40     mumps_.sym = 2;//general symmetric matrix
41     mumps_.comm_fortran = USE_COMM_WORLD;
42     int myid, ierr;
43     int justName;
44     ierr = MPI_Init(&justName, NULL);
45     ierr = MPI_Comm_rank(MPI_COMM_WORLD, &myid);
46     assert (!ierr);
47     dmumps_c(&mumps_);
48#define ICNTL(I) icntl[(I)-1] /* macro s.t. indices match documentation */
49#define CNTL(I) cntl[(I)-1] /* macro s.t. indices match documentation */
50     mumps_.ICNTL(5) = 1; // say compressed format
51     mumps_.ICNTL(4) = 2; // log messages
52     mumps_.ICNTL(24) = 1; // Deal with zeros on diagonal
53     mumps_.CNTL(3) = 1.0e-20; // drop if diagonal less than this
54     // output off
55     mumps_.ICNTL(1) = -1;
56     mumps_.ICNTL(2) = -1;
57     mumps_.ICNTL(3) = -1;
58     mumps_.ICNTL(4) = 0;
59}
60
61//-------------------------------------------------------------------
62// Copy constructor
63//-------------------------------------------------------------------
64ClpCholeskyMumps::ClpCholeskyMumps (const ClpCholeskyMumps & rhs)
65     : ClpCholeskyBase(rhs)
66{
67     abort();
68}
69
70
71//-------------------------------------------------------------------
72// Destructor
73//-------------------------------------------------------------------
74ClpCholeskyMumps::~ClpCholeskyMumps ()
75{
76     mumps_.job = JOB_END;
77     dmumps_c(&mumps_); /* Terminate instance */
78     MPI_Finalize();
79}
80
81//----------------------------------------------------------------
82// Assignment operator
83//-------------------------------------------------------------------
84ClpCholeskyMumps &
85ClpCholeskyMumps::operator=(const ClpCholeskyMumps& rhs)
86{
87     if (this != &rhs) {
88          ClpCholeskyBase::operator=(rhs);
89          abort();
90     }
91     return *this;
92}
93//-------------------------------------------------------------------
94// Clone
95//-------------------------------------------------------------------
96ClpCholeskyBase * ClpCholeskyMumps::clone() const
97{
98     return new ClpCholeskyMumps(*this);
99}
100/* Orders rows and saves pointer to matrix.and model */
101int
102ClpCholeskyMumps::order(ClpInterior * model)
103{
104     numberRows_ = model->numberRows();
105     if (doKKT_) {
106          numberRows_ += numberRows_ + model->numberColumns();
107          printf("finish coding MUMPS KKT!\n");
108          abort();
109     }
110     rowsDropped_ = new char [numberRows_];
111     memset(rowsDropped_, 0, numberRows_);
112     numberRowsDropped_ = 0;
113     model_ = model;
114     rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
115     const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
116     const int * columnLength = model_->clpMatrix()->getVectorLengths();
117     const int * row = model_->clpMatrix()->getIndices();
118     const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
119     const int * rowLength = rowCopy_->getVectorLengths();
120     const int * column = rowCopy_->getIndices();
121     // We need two arrays for counts
122     int * which = new int [numberRows_];
123     int * used = new int[numberRows_+1];
124     CoinZeroN(used, numberRows_);
125     int iRow;
126     sizeFactor_ = 0;
127     for (iRow = 0; iRow < numberRows_; iRow++) {
128          int number = 1;
129          // make sure diagonal exists
130          which[0] = iRow;
131          used[iRow] = 1;
132          if (!rowsDropped_[iRow]) {
133               CoinBigIndex startRow = rowStart[iRow];
134               CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
135               for (CoinBigIndex k = startRow; k < endRow; k++) {
136                    int iColumn = column[k];
137                    CoinBigIndex start = columnStart[iColumn];
138                    CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
139                    for (CoinBigIndex j = start; j < end; j++) {
140                         int jRow = row[j];
141                         if (jRow >= iRow && !rowsDropped_[jRow]) {
142                              if (!used[jRow]) {
143                                   used[jRow] = 1;
144                                   which[number++] = jRow;
145                              }
146                         }
147                    }
148               }
149               sizeFactor_ += number;
150               int j;
151               for (j = 0; j < number; j++)
152                    used[which[j]] = 0;
153          }
154     }
155     delete [] which;
156     // NOT COMPRESSED FOR NOW ??? - Space for starts
157     mumps_.ICNTL(5) = 0; // say NOT compressed format
158     try {
159          choleskyStart_ = new CoinBigIndex[numberRows_+1+sizeFactor_];
160     } catch (...) {
161          // no memory
162          return -1;
163     }
164     // Now we have size - create arrays and fill in
165     try {
166          choleskyRow_ = new int [sizeFactor_];
167     } catch (...) {
168          // no memory
169          delete [] choleskyStart_;
170          choleskyStart_ = NULL;
171          return -1;
172     }
173     try {
174          sparseFactor_ = new double[sizeFactor_];
175     } catch (...) {
176          // no memory
177          delete [] choleskyRow_;
178          choleskyRow_ = NULL;
179          delete [] choleskyStart_;
180          choleskyStart_ = NULL;
181          return -1;
182     }
183
184     sizeFactor_ = 0;
185     which = choleskyRow_;
186     for (iRow = 0; iRow < numberRows_; iRow++) {
187          int number = 1;
188          // make sure diagonal exists
189          which[0] = iRow;
190          used[iRow] = 1;
191          choleskyStart_[iRow] = sizeFactor_;
192          if (!rowsDropped_[iRow]) {
193               CoinBigIndex startRow = rowStart[iRow];
194               CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
195               for (CoinBigIndex k = startRow; k < endRow; k++) {
196                    int iColumn = column[k];
197                    CoinBigIndex start = columnStart[iColumn];
198                    CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
199                    for (CoinBigIndex j = start; j < end; j++) {
200                         int jRow = row[j];
201                         if (jRow >= iRow && !rowsDropped_[jRow]) {
202                              if (!used[jRow]) {
203                                   used[jRow] = 1;
204                                   which[number++] = jRow;
205                              }
206                         }
207                    }
208               }
209               sizeFactor_ += number;
210               int j;
211               for (j = 0; j < number; j++)
212                    used[which[j]] = 0;
213               // Sort
214               std::sort(which, which + number);
215               // move which on
216               which += number;
217          }
218     }
219     choleskyStart_[numberRows_] = sizeFactor_;
220     delete [] used;
221     permuteInverse_ = new int [numberRows_];
222     permute_ = new int[numberRows_];
223     // To Fortran and fake
224     for (iRow = 0; iRow < numberRows_ + 1; iRow++) {
225          int k = choleskyStart_[iRow];
226          int kEnd = choleskyStart_[iRow+1];
227          k += numberRows_ + 1;
228          kEnd += numberRows_ + 1;
229          for (; k < kEnd; k++)
230               choleskyStart_[k] = iRow + 1;
231          choleskyStart_[iRow]++;
232     }
233     mumps_.nz = sizeFactor_;
234     mumps_.irn = choleskyStart_ + numberRows_ + 1;
235     mumps_.jcn = choleskyRow_;
236     mumps_.a = NULL;
237     for (CoinBigIndex i = 0; i < sizeFactor_; i++) {
238          choleskyRow_[i]++;
239#ifndef NDEBUG
240          assert (mumps_.irn[i] >= 1 && mumps_.irn[i] <= numberRows_);
241          assert (mumps_.jcn[i] >= 1 && mumps_.jcn[i] <= numberRows_);
242#endif
243     }
244     // validate
245     //mumps code here
246     mumps_.n = numberRows_;
247     mumps_.nelt = numberRows_;
248     mumps_.eltptr = choleskyStart_;
249     mumps_.eltvar = choleskyRow_;
250     mumps_.a_elt = NULL;
251     mumps_.rhs = NULL;
252     mumps_.job = 1; // order
253     dmumps_c(&mumps_);
254     mumps_.a = sparseFactor_;
255     if (mumps_.infog[0]) {
256          printf("MUMPS ordering failed -error %d %d\n",
257                 mumps_.infog[0], mumps_.infog[1]);
258          return 1;
259     } else {
260          double size = mumps_.infog[19];
261          if (size < 0)
262               size *= -1000000;
263          printf("%g nonzeros, flop count %g\n", size, mumps_.rinfog[0]);
264     }
265     for (iRow = 0; iRow < numberRows_; iRow++) {
266          permuteInverse_[iRow] = iRow;
267          permute_[iRow] = iRow;
268     }
269     return 0;
270}
271/* Does Symbolic factorization given permutation.
272   This is called immediately after order.  If user provides this then
273   user must provide factorize and solve.  Otherwise the default factorization is used
274   returns non-zero if not enough memory */
275int
276ClpCholeskyMumps::symbolic()
277{
278     return 0;
279}
280/* Factorize - filling in rowsDropped and returning number dropped */
281int
282ClpCholeskyMumps::factorize(const double * diagonal, int * rowsDropped)
283{
284     const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
285     const int * columnLength = model_->clpMatrix()->getVectorLengths();
286     const int * row = model_->clpMatrix()->getIndices();
287     const double * element = model_->clpMatrix()->getElements();
288     const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
289     const int * rowLength = rowCopy_->getVectorLengths();
290     const int * column = rowCopy_->getIndices();
291     const double * elementByRow = rowCopy_->getElements();
292     int numberColumns = model_->clpMatrix()->getNumCols();
293     int iRow;
294     double * work = new double[numberRows_];
295     CoinZeroN(work, numberRows_);
296     const double * diagonalSlack = diagonal + numberColumns;
297     int newDropped = 0;
298     double largest;
299     //double smallest;
300     //perturbation
301     double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
302     perturbation = 0.0;
303     perturbation = perturbation * perturbation;
304     if (perturbation > 1.0) {
305#ifdef COIN_DEVELOP
306          //if (model_->model()->logLevel()&4)
307          std::cout << "large perturbation " << perturbation << std::endl;
308#endif
309          perturbation = sqrt(perturbation);;
310          perturbation = 1.0;
311     }
312     double delta2 = model_->delta(); // add delta*delta to diagonal
313     delta2 *= delta2;
314     for (iRow = 0; iRow < numberRows_; iRow++) {
315          double * put = sparseFactor_ + choleskyStart_[iRow] - 1; // Fortran
316          int * which = choleskyRow_ + choleskyStart_[iRow] - 1; // Fortran
317          int number = choleskyStart_[iRow+1] - choleskyStart_[iRow];
318          if (!rowLength[iRow])
319               rowsDropped_[iRow] = 1;
320          if (!rowsDropped_[iRow]) {
321               CoinBigIndex startRow = rowStart[iRow];
322               CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
323               work[iRow] = diagonalSlack[iRow] + delta2;
324               for (CoinBigIndex k = startRow; k < endRow; k++) {
325                    int iColumn = column[k];
326                    if (!whichDense_ || !whichDense_[iColumn]) {
327                         CoinBigIndex start = columnStart[iColumn];
328                         CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
329                         double multiplier = diagonal[iColumn] * elementByRow[k];
330                         for (CoinBigIndex j = start; j < end; j++) {
331                              int jRow = row[j];
332                              if (jRow >= iRow && !rowsDropped_[jRow]) {
333                                   double value = element[j] * multiplier;
334                                   work[jRow] += value;
335                              }
336                         }
337                    }
338               }
339               int j;
340               for (j = 0; j < number; j++) {
341                    int jRow = which[j] - 1; // to Fortran
342                    put[j] = work[jRow];
343                    work[jRow] = 0.0;
344               }
345          } else {
346               // dropped
347               int j;
348               for (j = 1; j < number; j++) {
349                    put[j] = 0.0;
350               }
351               put[0] = 1.0;
352          }
353     }
354     //check sizes
355     double largest2 = maximumAbsElement(sparseFactor_, sizeFactor_);
356     largest2 *= 1.0e-20;
357     largest = CoinMin(largest2, 1.0e-11);
358     int numberDroppedBefore = 0;
359     for (iRow = 0; iRow < numberRows_; iRow++) {
360          int dropped = rowsDropped_[iRow];
361          // Move to int array
362          rowsDropped[iRow] = dropped;
363          if (!dropped) {
364               CoinBigIndex start = choleskyStart_[iRow] - 1; // to Fortran
365               double diagonal = sparseFactor_[start];
366               if (diagonal > largest2) {
367                    sparseFactor_[start] = CoinMax(diagonal, 1.0e-10);
368               } else {
369                    sparseFactor_[start] = CoinMax(diagonal, 1.0e-10);
370                    rowsDropped[iRow] = 2;
371                    numberDroppedBefore++;
372               }
373          }
374     }
375     delete [] work;
376     // code here
377     mumps_.a_elt = sparseFactor_;
378     mumps_.rhs = NULL;
379     mumps_.job = 2; // factorize
380     dmumps_c(&mumps_);
381     if (mumps_.infog[0]) {
382          printf("MUMPS factorization failed -error %d %d\n",
383                 mumps_.infog[0], mumps_.infog[1]);
384     }
385     choleskyCondition_ = 1.0;
386     bool cleanCholesky;
387     if (model_->numberIterations() < 2000)
388          cleanCholesky = true;
389     else
390          cleanCholesky = false;
391     if (cleanCholesky) {
392          //drop fresh makes some formADAT easier
393          //int oldDropped=numberRowsDropped_;
394          if (newDropped || numberRowsDropped_) {
395               //std::cout <<"Rank "<<numberRows_-newDropped<<" ( "<<
396               //  newDropped<<" dropped)";
397               //if (newDropped>oldDropped)
398               //std::cout<<" ( "<<newDropped-oldDropped<<" dropped this time)";
399               //std::cout<<std::endl;
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] = 0;
408                    }
409               }
410               numberRowsDropped_ = newDropped;
411               newDropped = -(2 + newDropped);
412          }
413     } else {
414          if (newDropped) {
415               newDropped = 0;
416               for (int i = 0; i < numberRows_; i++) {
417                    char dropped = rowsDropped[i];
418                    rowsDropped_[i] = dropped;
419                    if (dropped == 2) {
420                         //dropped this time
421                         rowsDropped[newDropped++] = i;
422                         rowsDropped_[i] = 1;
423                    }
424               }
425          }
426          numberRowsDropped_ += newDropped;
427          if (numberRowsDropped_ && 0) {
428               std::cout << "Rank " << numberRows_ - numberRowsDropped_ << " ( " <<
429                         numberRowsDropped_ << " dropped)";
430               if (newDropped) {
431                    std::cout << " ( " << newDropped << " dropped this time)";
432               }
433               std::cout << std::endl;
434          }
435     }
436     status_ = 0;
437     return newDropped;
438}
439/* Uses factorization to solve. */
440void
441ClpCholeskyMumps::solve (double * region)
442{
443     mumps_.rhs = region;
444     mumps_.job = 3; // solve
445     dmumps_c(&mumps_);
446}
447#endif
Note: See TracBrowser for help on using the repository browser.