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

Last change on this file since 2030 was 1929, checked in by stefan, 6 years ago

fix compiler (gcc 4.6.2) warnings in optimized mode, mainly about unused variables

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