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

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