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

Last change on this file since 553 was 553, checked in by andreasw, 14 years ago

added latest changes proposed by Olaf Schenk in PardisoInterface?

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