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

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

out some printf statements

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