source: branches/dev/Algorithm/LinearSolvers/IpPardisoSolverInterface.cpp @ 547

Last change on this file since 547 was 547, checked in by andreasw, 14 years ago
  • corrected infinite-loop bug in monotone mu mupdate
  • perturbation for constraints now as in paper (with exponent for mu)
  • some changes in Pardiso interface
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.1 KB
Line 
1// Copyright (C) 2005 International Business Machines and others.
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// $Id: IpPardisoSolverInterface.cpp 547 2005-10-24 15:24:19Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2005-03-17
8//
9//           Olaf Schenk                      Univ of Basel 2005-09-20
10//                  - changed options, added PHASE_ flag
11
12#include "IpPardisoSolverInterface.hpp"
13
14
15#ifdef HAVE_CMATH
16# include <cmath>
17#else
18# ifdef HAVE_MATH_H
19#  include <math.h>
20# else
21#  error "don't have header file for math"
22# endif
23#endif
24
25/** Prototypes for Pardiso's subroutines */
26extern "C"
27{
28  void F77_FUNC(pardisoinit,PARDISOINIT)(void* PT, const ipfint* MTYPE,
29                                         ipfint* IPARM);
30  void F77_FUNC(pardiso,PARDISO)(void** PT, const ipfint* MAXFCT,
31                                 const ipfint* MNUM, const ipfint* MTYPE,
32                                 const ipfint* PHASE, const ipfint* N,
33                                 const double* A, const ipfint* IA,
34                                 const ipfint* JA, const ipfint* PERM,
35                                 const ipfint* NRHS, ipfint* IPARM,
36                                 const ipfint* MSGLVL, double* B, double* X,
37                                 ipfint* ERROR);
38}
39
40namespace Ipopt
41{
42#ifdef IP_DEBUG
43  static const Index dbg_verbosity = 0;
44#endif
45
46  PardisoSolverInterface::PardisoSolverInterface()
47      :
48      a_(NULL),
49      initialized_(false),
50
51      MAXFCT_(1),
52      MNUM_(1),
53      MTYPE_(-2),
54      MSGLVL_(0)
55  {
56    DBG_START_METH("PardisoSolverInterface::PardisoSolverInterface()",dbg_verbosity);
57
58    PT_ = new void*[64];
59    IPARM_ = new ipfint[64];
60  }
61
62  PardisoSolverInterface::~PardisoSolverInterface()
63  {
64    DBG_START_METH("PardisoSolverInterface::~PardisoSolverInterface()",
65                   dbg_verbosity);
66
67    // Tell Pardiso to release all memory
68    if (initialized_) {
69      ipfint PHASE = -1;
70      ipfint N = dim_;
71      ipfint NRHS = 0;
72      ipfint ERROR;
73      F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_, &PHASE, &N,
74                                NULL, NULL, NULL, NULL, &NRHS, IPARM_,
75                                &MSGLVL_, NULL, NULL, &ERROR);
76      DBG_ASSERT(ERROR==0);
77    }
78
79    delete[] PT_;
80    delete[] IPARM_;
81    delete[] a_;
82  }
83
84  bool PardisoSolverInterface::InitializeImpl(const OptionsList& options,
85      const std::string& prefix)
86  {
87    // Number value = 0.0;
88
89    // Tell Pardiso to release all memory if it had been used before
90    if (initialized_) {
91      ipfint PHASE = -1;
92      ipfint N = dim_;
93      ipfint NRHS = 0;
94      ipfint ERROR;
95      F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_, &PHASE, &N,
96                                NULL, NULL, NULL, NULL, &NRHS, IPARM_,
97                                &MSGLVL_, NULL, NULL, &ERROR);
98      DBG_ASSERT(ERROR==0);
99    }
100
101    // Reset all private data
102    dim_=0;
103    nonzeros_=0;
104    have_symbolic_factorization_=false;
105    initialized_=false;
106    delete[] a_;
107
108    // Call Pardiso's initialization routine
109    IPARM_[0] = 0;  // Tell it to fill IPARM with default values(?)
110    F77_FUNC(pardisoinit,PARDISOINIT)(PT_, &MTYPE_, IPARM_);
111
112    // Set some parameters for Pardiso
113    IPARM_[0] = 1;  // Don't use the default values
114    IPARM_[2] = 1;  // Only one CPU for now
115    IPARM_[5] = 1;  // Overwrite right-hand side
116    // ToDo: decide if we need iterative refinement in Pardiso.  For
117    // now, switch it off ?
118    IPARM_[7] = 0;
119
120    // Options suggested by Olaf Schenk
121    IPARM_[9] = 12;
122    IPARM_[10] = 1;
123    // Matching information:  IPARM_[12] = 1 seems ok, but results in a
124    // large number of pivot perturbation
125    // Matching information:  IPARM_[12] = 2 robust,  but more  expensive method
126    IPARM_[12] = 2;
127
128    IPARM_[20] = 1;
129
130    return true;
131  }
132
133  ESymSolverStatus PardisoSolverInterface::MultiSolve(bool new_matrix,
134      const Index* ia,
135      const Index* ja,
136      Index nrhs,
137      double* rhs_vals,
138      bool check_NegEVals,
139      Index numberOfNegEVals)
140  {
141    DBG_START_METH("PardisoSolverInterface::MultiSolve",dbg_verbosity);
142    DBG_ASSERT(!check_NegEVals || ProvidesInertia());
143    DBG_ASSERT(initialized_);
144
145    // check if a factorization has to be done
146    if (new_matrix) {
147      // perform the factorization
148      ESymSolverStatus retval;
149      retval = Factorization(ia, ja, check_NegEVals, numberOfNegEVals);
150      if (retval!=SYMSOLVER_SUCCESS) {
151        DBG_PRINT((1, "FACTORIZATION FAILED!\n"));
152        return retval;  // Matrix singular or error occurred
153      }
154    }
155
156    // do the solve
157    return Solve(ia, ja, nrhs, rhs_vals);
158  }
159
160  double* PardisoSolverInterface::GetValuesArrayPtr()
161  {
162    DBG_ASSERT(initialized_);
163    DBG_ASSERT(a_);
164    return a_;
165  }
166
167  /** Initialize the local copy of the positions of the nonzero
168      elements */
169  ESymSolverStatus PardisoSolverInterface::InitializeStructure
170  (Index dim, Index nonzeros,
171   const Index* ia,
172   const Index* ja)
173  {
174    DBG_START_METH("PardisoSolverInterface::InitializeStructure",dbg_verbosity);
175    dim_ = dim;
176    nonzeros_ = nonzeros;
177
178    // Make space for storing the matrix elements
179    delete[] a_;
180    a_ = new double[nonzeros_];
181
182    // Do the symbolic facotrization
183    ESymSolverStatus retval = SymbolicFactorization(ia, ja);
184    if (retval != SYMSOLVER_SUCCESS) {
185      return retval;
186    }
187
188    initialized_ = true;
189
190    return retval;
191  }
192
193  ESymSolverStatus
194  PardisoSolverInterface::SymbolicFactorization(const Index* ia,
195      const Index* ja)
196  {
197    DBG_START_METH("PardisoSolverInterface::SymbolicFactorization",
198                   dbg_verbosity);
199
200    // Since Pardiso requires the values of the nonzeros of the matrix
201    // for an efficient symbolic factorization, we postpone that task
202    // until the first call of Factorize.  All we do here is to reset
203    // the flag (in case this interface is called for a matrix with a
204    // new structure).
205
206    have_symbolic_factorization_ = false;
207
208    return SYMSOLVER_SUCCESS;
209  }
210
211  ESymSolverStatus
212  PardisoSolverInterface::Factorization(const Index* ia,
213                                        const Index* ja,
214                                        bool check_NegEVals,
215                                        Index numberOfNegEVals)
216  {
217    DBG_START_METH("PardisoSolverInterface::Factorization",dbg_verbosity);
218
219    // Call Pardiso to do the factorization
220    ipfint N = dim_;
221    ipfint PHASE;
222    ipfint PERM;   // This should not be accessed by Pardiso
223    ipfint NRHS = 0;
224    double B;  // This should not be accessed by Pardiso in factorization
225    // phase
226    double X;  // This should not be accessed by Pardiso in factorization
227    // phase
228    ipfint ERROR;
229
230    if (getenv ("IPOPT_WRITE_MAT")) {
231
232      // Dump matrix to file...
233      ipfint  NNZ = ia[N]-1;
234      ipfint   i;
235
236      if (IpData().iter_count() != debug_last_iter_) {
237        debug_cnt_ = 0;
238      }
239
240      /* Write header */
241      FILE    *mat_file;
242      char     mat_name[128];
243      char     mat_pref[32];
244
245      if (getenv ("IPOPT_WRITE_PREFIX"))
246        strcpy (mat_pref, getenv ("IPOPT_WRITE_PREFIX"));
247      else
248        strcpy (mat_pref, "mat-ipopt");
249
250      sprintf (mat_name, "%s_%03d-%02d.iajaa", mat_pref, IpData().iter_count(), debug_cnt_);
251      mat_file = fopen (mat_name, "w");
252
253      fprintf (mat_file, "%d\n", N);
254      fprintf (mat_file, "%d\n", NNZ);
255
256      /* Write ia's */
257      for (i = 0; i < N+1; i++)
258        fprintf (mat_file, "%d\n", ia[i]);
259
260      /* Write ja's */
261      for (i = 0; i < NNZ; i++)
262        fprintf (mat_file, "%d\n", ja[i]);
263
264      /* Write values */
265      for (i = 0; i < NNZ; i++)
266        fprintf (mat_file, "%32.24e\n", a_[i]);
267
268      // /* Write RHS */
269      // for (i = 0; i < N; i++)
270      //        fprintf (mat_file, "%32.24e\n", rhs_vals[i]);
271
272      fclose (mat_file);
273
274      debug_last_iter_ = IpData().iter_count();
275      debug_cnt_++;
276    }
277
278    bool done = false;
279    bool just_performed_symbolic_factorization = false;
280
281    while (!done) {
282      if (!have_symbolic_factorization_) {
283        IpData().TimingStats().LinearSystemSymbolicFactorization().Start();
284        PHASE = 11;
285        Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
286                       "Calling Pardiso for symbolic factorization.\n");
287        F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_,
288                                  &PHASE, &N, a_, ia, ja, &PERM,
289                                  &NRHS, IPARM_, &MSGLVL_, &B, &X,
290                                  &ERROR);
291        IpData().TimingStats().LinearSystemSymbolicFactorization().End();
292        if (ERROR!=0 ) {
293          Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
294                         "Error in Pardiso during symbolic factorization phase.  ERROR = %d.\n", ERROR);
295          return SYMSOLVER_FATAL_ERROR;
296        }
297        have_symbolic_factorization_ = true;
298        just_performed_symbolic_factorization = true;
299      }
300
301      PHASE = 22;
302
303      IpData().TimingStats().LinearSystemFactorization().Start();
304      Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
305                     "Calling Pardiso for factorization.\n");
306      F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_,
307                                &PHASE, &N, a_, ia, ja, &PERM,
308                                &NRHS, IPARM_, &MSGLVL_, &B, &X,
309                                &ERROR);
310      IpData().TimingStats().LinearSystemFactorization().End();
311
312      if (ERROR==-4) {
313        // I think this means that the matrix is singular
314        // OLAF said that this will never happen (ToDo)
315        return SYMSOLVER_SINGULAR;
316      }
317      else if (ERROR!=0 ) {
318        Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
319                       "Error in Pardiso during factorization phase.  ERROR = %d.\n", ERROR);
320        return SYMSOLVER_FATAL_ERROR;
321      }
322
323      negevals_ = IPARM_[22];
324      if (IPARM_[13] != 0) {
325        Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
326                       "Number of perturbed pivots in factorization phase = %d.\n", IPARM_[13]);
327        // trigger a new symblic factorization for the next factorize
328        // call, even if the current inertia estimate is correct
329        // OLAF MICHAEL : Maybe we need to discuss this
330        have_symbolic_factorization_ = false;
331        // We assume now that if there was just a symbolic
332        // factorization and we still have perturbed pivots, that the
333        // system is actually singular
334        if (just_performed_symbolic_factorization) {
335          IpData().Append_info_string("Ps");
336          //break; // Declaring this system singular doesn't seem to be
337                 // working well???
338          return SYMSOLVER_SINGULAR;
339          }
340        IpData().Append_info_string("Pp");
341        /*
342        // If we there were perturbed pivots, and the inertia of the
343        // system is correct, and we haven't redone the symblic
344        // factorization during this call yet, trigger a new
345        // factorization after a symblic factorization
346        done = (just_performed_symbolic_factorization || !check_NegEVals ||
347                numberOfNegEVals==negevals_);
348        */
349        done = false;
350      }
351      else {
352        done = true;
353      }
354    }
355
356    DBG_ASSERT(IPARM_[21]+IPARM_[22] == dim_);
357
358    // Check whether the number of negative eigenvalues matches the requested
359    // count
360    if (check_NegEVals && (numberOfNegEVals!=negevals_)) {
361      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
362                     "Wrong inertia: required are %d, but we got %d.\n",
363                     numberOfNegEVals, negevals_);
364      return SYMSOLVER_WRONG_INERTIA;
365    }
366
367    return SYMSOLVER_SUCCESS;
368  }
369
370  ESymSolverStatus PardisoSolverInterface::Solve(const Index* ia,
371      const Index* ja,
372      Index nrhs,
373      double *rhs_vals)
374  {
375    DBG_START_METH("PardisoSolverInterface::Solve",dbg_verbosity);
376
377    IpData().TimingStats().LinearSystemBackSolve().Start();
378    // Call Pardiso to do the solve for the given right-hand sides
379    ipfint PHASE = 33;
380    ipfint N = dim_;
381    ipfint PERM;   // This should not be accessed by Pardiso
382    ipfint NRHS = nrhs;
383    double* X = new double[nrhs*dim_];
384    ipfint ERROR;
385
386    F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_,
387                              &PHASE, &N, a_, ia, ja, &PERM,
388                              &NRHS, IPARM_, &MSGLVL_, rhs_vals, X,
389                              &ERROR);
390
391    // The following delete was missing (memory leak?)
392    delete [] X; /* OLAF/MICHAEL: do we really need X? */
393
394    if (IPARM_[6] != 0) {
395      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
396                     "Number of iterative refinement steps = %d.\n", IPARM_[6]);
397      IpData().Append_info_string("Pi");
398    }
399
400    IpData().TimingStats().LinearSystemBackSolve().End();
401    if (ERROR!=0 ) {
402      Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
403                     "Error in Pardiso during solve phase.  ERROR = %d.\n", ERROR);
404      return SYMSOLVER_FATAL_ERROR;
405    }
406    return SYMSOLVER_SUCCESS;
407  }
408
409  Index PardisoSolverInterface::NumberOfNegEVals() const
410  {
411    DBG_START_METH("PardisoSolverInterface::NumberOfNegEVals",dbg_verbosity);
412    DBG_ASSERT(negevals_>=0);
413    return negevals_;
414  }
415
416  bool PardisoSolverInterface::IncreaseQuality()
417  {
418    // At the moment, I don't see how we could tell Pardiso to do better
419    // (maybe switch from IPARM[20]=1 to IPARM[20]=2?)
420    return false;
421  }
422
423} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.