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

Last change on this file since 2470 was 2385, checked in by unxusr, 9 months ago

formatting

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