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

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

astyle

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.2 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 534 2005-10-04 15:48:18Z 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    IPARM_[12] = 1;
124
125    IPARM_[20] = 1;
126
127    return true;
128  }
129
130  ESymSolverStatus PardisoSolverInterface::MultiSolve(bool new_matrix,
131      const Index* ia,
132      const Index* ja,
133      Index nrhs,
134      double* rhs_vals,
135      bool check_NegEVals,
136      Index numberOfNegEVals)
137  {
138    DBG_START_METH("PardisoSolverInterface::MultiSolve",dbg_verbosity);
139    DBG_ASSERT(!check_NegEVals || ProvidesInertia());
140    DBG_ASSERT(initialized_);
141
142    // check if a factorization has to be done
143    if (new_matrix) {
144      // perform the factorization
145      ESymSolverStatus retval;
146      retval = Factorization(ia, ja, check_NegEVals, numberOfNegEVals);
147      if (retval!=SYMSOLVER_SUCCESS) {
148        DBG_PRINT((1, "FACTORIZATION FAILED!\n"));
149        return retval;  // Matrix singular or error occurred
150      }
151    }
152
153    // do the solve
154    return Solve(ia, ja, nrhs, rhs_vals);
155  }
156
157  double* PardisoSolverInterface::GetValuesArrayPtr()
158  {
159    DBG_ASSERT(initialized_);
160    DBG_ASSERT(a_);
161    return a_;
162  }
163
164  /** Initialize the local copy of the positions of the nonzero
165      elements */
166  ESymSolverStatus PardisoSolverInterface::InitializeStructure
167  (Index dim, Index nonzeros,
168   const Index* ia,
169   const Index* ja)
170  {
171    DBG_START_METH("PardisoSolverInterface::InitializeStructure",dbg_verbosity);
172    dim_ = dim;
173    nonzeros_ = nonzeros;
174
175    // Make space for storing the matrix elements
176    delete[] a_;
177    a_ = new double[nonzeros_];
178
179    // Do the symbolic facotrization
180    ESymSolverStatus retval = SymbolicFactorization(ia, ja);
181    if (retval != SYMSOLVER_SUCCESS) {
182      return retval;
183    }
184
185    initialized_ = true;
186
187    return retval;
188  }
189
190  ESymSolverStatus
191  PardisoSolverInterface::SymbolicFactorization(const Index* ia,
192      const Index* ja)
193  {
194    DBG_START_METH("PardisoSolverInterface::SymbolicFactorization",
195                   dbg_verbosity);
196
197    // Since Pardiso requires the values of the nonzeros of the matrix
198    // for an efficient symbolic factorization, we postpone that task
199    // until the first call of Factorize.  All we do here is to reset
200    // the flag (in case this interface is called for a matrix with a
201    // new structure).
202
203    have_symbolic_factorization_ = false;
204
205    return SYMSOLVER_SUCCESS;
206  }
207
208  ESymSolverStatus
209  PardisoSolverInterface::Factorization(const Index* ia,
210                                        const Index* ja,
211                                        bool check_NegEVals,
212                                        Index numberOfNegEVals)
213  {
214    DBG_START_METH("PardisoSolverInterface::Factorization",dbg_verbosity);
215
216    // Call Pardiso to do the factorization
217    ipfint N = dim_;
218    ipfint PHASE;
219    ipfint PERM;   // This should not be accessed by Pardiso
220    ipfint NRHS = 0;
221    double B;  // This should not be accessed by Pardiso in factorization
222    // phase
223    double X;  // This should not be accessed by Pardiso in factorization
224    // phase
225    ipfint ERROR;
226
227    if (getenv ("IPOPT_WRITE_MAT")) {
228
229      // Dump matrix to file...
230      ipfint  NNZ = ia[N]-1;
231      ipfint   i;
232
233      if (IpData().iter_count() != debug_last_iter_) {
234        debug_cnt_ = 0;
235      }
236
237      /* Write header */
238      FILE    *mat_file;
239      char     mat_name[128];
240      char     mat_pref[32];
241
242      if (getenv ("IPOPT_WRITE_PREFIX"))
243        strcpy (mat_pref, getenv ("IPOPT_WRITE_PREFIX"));
244      else
245        strcpy (mat_pref, "mat-ipopt");
246
247      sprintf (mat_name, "%s_%03d-%02d.iajaa", mat_pref, IpData().iter_count(), debug_cnt_);
248      mat_file = fopen (mat_name, "w");
249
250      fprintf (mat_file, "%d\n", N);
251      fprintf (mat_file, "%d\n", NNZ);
252
253      /* Write ia's */
254      for (i = 0; i < N+1; i++)
255        fprintf (mat_file, "%d\n", ia[i]);
256
257      /* Write ja's */
258      for (i = 0; i < NNZ; i++)
259        fprintf (mat_file, "%d\n", ja[i]);
260
261      /* Write values */
262      for (i = 0; i < NNZ; i++)
263        fprintf (mat_file, "%32.24e\n", a_[i]);
264
265      // /* Write RHS */
266      // for (i = 0; i < N; i++)
267      //        fprintf (mat_file, "%32.24e\n", rhs_vals[i]);
268
269      fclose (mat_file);
270
271      debug_last_iter_ = IpData().iter_count();
272      debug_cnt_++;
273    }
274
275    bool done = false;
276    bool just_performed_symbolic_factorization = false;
277
278    while (!done) {
279      if (!have_symbolic_factorization_) {
280        IpData().TimingStats().LinearSystemSymbolicFactorization().Start();
281        PHASE = 11;
282        F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_,
283                                  &PHASE, &N, a_, ia, ja, &PERM,
284                                  &NRHS, IPARM_, &MSGLVL_, &B, &X,
285                                  &ERROR);
286        IpData().TimingStats().LinearSystemSymbolicFactorization().End();
287        if (ERROR!=0 ) {
288          Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
289                         "Error in Pardiso during symbolic factorization phase.  ERROR = %d.\n", ERROR);
290          return SYMSOLVER_FATAL_ERROR;
291        }
292        have_symbolic_factorization_ = true;
293        just_performed_symbolic_factorization = true;
294      }
295
296      PHASE = 22;
297
298      IpData().TimingStats().LinearSystemFactorization().Start();
299      F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_,
300                                &PHASE, &N, a_, ia, ja, &PERM,
301                                &NRHS, IPARM_, &MSGLVL_, &B, &X,
302                                &ERROR);
303      IpData().TimingStats().LinearSystemFactorization().End();
304
305      if (ERROR==-4) {
306        // I think this means that the matrix is singular
307        // OLAF said that this will never happen (ToDo)
308        return SYMSOLVER_SINGULAR;
309      }
310      else if (ERROR!=0 ) {
311        Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
312                       "Error in Pardiso during factorization phase.  ERROR = %d.\n", ERROR);
313        return SYMSOLVER_FATAL_ERROR;
314      }
315
316      negevals_ = IPARM_[22];
317      if (IPARM_[13] != 0) {
318        Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
319                       "Number of perturbed pivots in factorization phase = %d.\n", IPARM_[13]);
320        IpData().Append_info_string("Pp");
321        // trigger a new symblic factorization for the next factorize
322        // call, even if the current inertia estimate is correct
323        // OLAF MICHAEL : Maybe we need to discuss this
324        have_symbolic_factorization_ = false;
325        // If we there were perturbed pivots, and the inertia of the
326        // system is correct, and we haven't redone the symblic
327        // factorization during this call yet, trigger a new
328        // factorization after a symblic factorization
329        done = (just_performed_symbolic_factorization || !check_NegEVals ||
330                numberOfNegEVals==negevals_);
331      }
332      else {
333        done = true;
334      }
335    }
336
337    DBG_ASSERT(IPARM_[21]+IPARM_[22] == dim_);
338
339    // Check whether the number of negative eigenvalues matches the requested
340    // count
341    if (check_NegEVals && (numberOfNegEVals!=negevals_)) {
342      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
343                     "Wrong inertia: required are %d, but we got %d.\n",
344                     numberOfNegEVals, negevals_);
345      return SYMSOLVER_WRONG_INERTIA;
346    }
347
348    return SYMSOLVER_SUCCESS;
349  }
350
351  ESymSolverStatus PardisoSolverInterface::Solve(const Index* ia,
352      const Index* ja,
353      Index nrhs,
354      double *rhs_vals)
355  {
356    DBG_START_METH("PardisoSolverInterface::Solve",dbg_verbosity);
357
358    IpData().TimingStats().LinearSystemBackSolve().Start();
359    // Call Pardiso to do the solve for the given right-hand sides
360    ipfint PHASE = 33;
361    ipfint N = dim_;
362    ipfint PERM;   // This should not be accessed by Pardiso
363    ipfint NRHS = nrhs;
364    double* X = new double[nrhs*dim_];
365    ipfint ERROR;
366
367    F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_,
368                              &PHASE, &N, a_, ia, ja, &PERM,
369                              &NRHS, IPARM_, &MSGLVL_, rhs_vals, X,
370                              &ERROR);
371
372    if (IPARM_[6] != 0) {
373      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
374                     "Number of iterative refinement steps = %d.\n", IPARM_[6]);
375      IpData().Append_info_string("Pi");
376    }
377
378    IpData().TimingStats().LinearSystemBackSolve().End();
379    if (ERROR!=0 ) {
380      Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
381                     "Error in Pardiso during solve phase.  ERROR = %d.\n", ERROR);
382      return SYMSOLVER_FATAL_ERROR;
383    }
384    return SYMSOLVER_SUCCESS;
385  }
386
387  Index PardisoSolverInterface::NumberOfNegEVals() const
388  {
389    DBG_START_METH("PardisoSolverInterface::NumberOfNegEVals",dbg_verbosity);
390    DBG_ASSERT(negevals_>=0);
391    return negevals_;
392  }
393
394  bool PardisoSolverInterface::IncreaseQuality()
395  {
396    // At the moment, I don't see how we could tell Pardiso to do better
397    // (maybe switch from IPARM[20]=1 to IPARM[20]=2?)
398    return false;
399  }
400
401} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.