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

Last change on this file since 2247 was 2247, checked in by forrest, 3 years ago

safer to use mpi.h not mumps_mpi.h

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