source: branches/dev/PDSystemImpl/PDFullSpace/IpPardisoSolverInterface.cpp @ 427

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

In case of wrong inertia: in the number of negative eigenvalues is TOO
SMALL, we now first try to fix that by increaseing pivot tolerance,
and otherwise pretend that the system is singular.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.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 427 2005-08-05 20:34:00Z andreasw $
6//
7// Authors:  Andreas Waechter     IBM    2005-03-17
8
9#include "IpPardisoSolverInterface.hpp"
10
11#ifdef OLD_C_HEADERS
12# include <math.h>
13#else
14# include <cmath>
15#endif
16
17/** Prototypes for Pardiso's subroutines */
18extern "C"
19{
20  void F77_FUNC(pardisoinit,PARDISOINIT)(void* PT, const ipfint* MTYPE,
21                                         ipfint* IPARM);
22  void F77_FUNC(pardiso,PARDISO)(void** PT, const ipfint* MAXFCT,
23                                 const ipfint* MNUM, const ipfint* MTYPE,
24                                 const ipfint* PHASE, const ipfint* N,
25                                 const double* A, const ipfint* IA,
26                                 const ipfint* JA, const ipfint* PERM,
27                                 const ipfint* NRHS, ipfint* IPARM,
28                                 const ipfint* MSGLVL, double* B, double* X,
29                                 ipfint* ERROR);
30}
31
32namespace Ipopt
33{
34
35  DBG_SET_VERBOSITY(0);
36
37  PardisoSolverInterface::PardisoSolverInterface()
38      :
39      dim_(0),
40      nonzeros_(0),
41      initialized_(false),
42      negevals_(-1),
43
44      MAXFCT_(1),
45      MNUM_(1),
46      MTYPE_(-2),
47      MSGLVL_(0),
48
49      a_(NULL)
50  {
51    DBG_START_METH("PardisoSolverInterface::PardisoSolverInterface()",dbg_verbosity);
52
53    PT_ = new void*[64];
54    IPARM_ = new ipfint[64];
55  }
56
57  PardisoSolverInterface::~PardisoSolverInterface()
58  {
59    DBG_START_METH("PardisoSolverInterface::~PardisoSolverInterface()",
60                   dbg_verbosity);
61
62    // Tell Pardiso to release all memory
63    if (initialized_) {
64      ipfint PHASE = -1;
65      ipfint N = dim_;
66      ipfint NRHS = 0;
67      ipfint ERROR;
68      F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_, &PHASE, &N,
69                                NULL, NULL, NULL, NULL, &NRHS, IPARM_,
70                                &MSGLVL_, NULL, NULL, &ERROR);
71      DBG_ASSERT(ERROR==0);
72    }
73
74    delete[] PT_;
75    delete[] IPARM_;
76    delete[] a_;
77  }
78
79  bool PardisoSolverInterface::InitializeImpl(const OptionsList& options,
80      const std::string& prefix)
81  {
82    Number value = 0.0;
83
84    // Tell Pardiso to release all memory if it had been used before
85    if (initialized_) {
86      ipfint PHASE = -1;
87      ipfint N = dim_;
88      ipfint NRHS = 0;
89      ipfint ERROR;
90      F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_, &PHASE, &N,
91                                NULL, NULL, NULL, NULL, &NRHS, IPARM_,
92                                &MSGLVL_, NULL, NULL, &ERROR);
93      DBG_ASSERT(ERROR==0);
94    }
95
96    // Reset all private data
97    dim_=0;
98    nonzeros_=0;
99    initialized_=false;
100    delete[] a_;
101
102    // Call Pardiso's initialization routine
103    IPARM_[0] = 0;  // Tell it to fill IPARM with default values(?)
104    F77_FUNC(pardisoinit,PARDISOINIT)(PT_, &MTYPE_, IPARM_);
105
106    // Set some parameters for Pardiso
107    IPARM_[0] = 1;  // Don't use the default values
108    IPARM_[2] = 1;  // Only one CPU for now
109    IPARM_[5] = 1;  // Overwrite right-hand side
110
111    // ToDo: decide if we need iterative refinement in Pardiso.  For
112    // now, switch it off ?
113    IPARM_[7] = 0;
114
115    // IPARM_[20] = 2;
116
117    return true;
118  }
119
120  ESymSolverStatus PardisoSolverInterface::MultiSolve(bool new_matrix,
121      const Index* ia,
122      const Index* ja,
123      Index nrhs,
124      double* rhs_vals,
125      bool check_NegEVals,
126      Index numberOfNegEVals)
127  {
128    DBG_START_METH("PardisoSolverInterface::MultiSolve",dbg_verbosity);
129    DBG_ASSERT(!check_NegEVals || ProvidesInertia());
130    DBG_ASSERT(initialized_);
131
132    // check if a factorization has to be done
133    if (new_matrix) {
134      // perform the factorization
135      ESymSolverStatus retval;
136      retval = Factorization(ia, ja, check_NegEVals, numberOfNegEVals);
137      if (retval!=SYMSOLVER_SUCCESS) {
138        DBG_PRINT((1, "FACTORIZATION FAILED!\n"));
139        return retval;  // Matrix singular or error occurred
140      }
141    }
142
143    // do the solve
144    return Solve(ia, ja, nrhs, rhs_vals);
145  }
146
147  double* PardisoSolverInterface::GetValuesArrayPtr()
148  {
149    DBG_ASSERT(initialized_);
150    DBG_ASSERT(a_);
151    return a_;
152  }
153
154  /** Initialize the local copy of the positions of the nonzero
155      elements */
156  ESymSolverStatus PardisoSolverInterface::InitializeStructure
157  (Index dim, Index nonzeros,
158   const Index* ia,
159   const Index* ja)
160  {
161    DBG_START_METH("PardisoSolverInterface::InitializeStructure",dbg_verbosity);
162    dim_ = dim;
163    nonzeros_ = nonzeros;
164
165    // Make space for storing the matrix elements
166    delete[] a_;
167    a_ = new double[nonzeros_];
168
169    // Do the symbolic facotrization
170    ESymSolverStatus retval = SymbolicFactorization(ia, ja);
171    if (retval != SYMSOLVER_SUCCESS) {
172      return retval;
173    }
174
175    initialized_ = true;
176
177    return retval;
178  }
179
180  ESymSolverStatus
181  PardisoSolverInterface::SymbolicFactorization(const Index* ia,
182      const Index* ja)
183  {
184    DBG_START_METH("PardisoSolverInterface::SymbolicFactorization",
185                   dbg_verbosity);
186
187    // Call Pardiso to do the analysis phase
188    ipfint PHASE = 11;
189    ipfint N = dim_;
190    ipfint PERM;   // This should not be accessed by Pardiso
191    ipfint NRHS = 0;
192    double B;  // This should not be accessed by Pardiso in analysis
193    // phase
194    double X;  // This should not be accessed by Pardiso in analysis
195    // phase
196    ipfint ERROR;
197
198    F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_,
199                              &PHASE, &N, a_, ia, ja, &PERM,
200                              &NRHS, IPARM_, &MSGLVL_, &B, &X,
201                              &ERROR);
202    if (ERROR!=0) {
203      Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
204                     "Error in Pardiso during analysis phase.  ERROR = %d.\n",
205                     ERROR);
206      return SYMSOLVER_FATAL_ERROR;
207    }
208
209    return SYMSOLVER_SUCCESS;
210  }
211
212  ESymSolverStatus
213  PardisoSolverInterface::Factorization(const Index* ia,
214                                        const Index* ja,
215                                        bool check_NegEVals,
216                                        Index numberOfNegEVals)
217  {
218    DBG_START_METH("PardisoSolverInterface::Factorization",dbg_verbosity);
219
220    // Call Pardiso to do the factorization
221    ipfint PHASE = 22;
222    ipfint N = dim_;
223    ipfint PERM;   // This should not be accessed by Pardiso
224    ipfint NRHS = 0;
225    double B;  // This should not be accessed by Pardiso in factorization
226    // phase
227    double X;  // This should not be accessed by Pardiso in factorization
228    // phase
229    ipfint ERROR;
230
231    F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_,
232                              &PHASE, &N, a_, ia, ja, &PERM,
233                              &NRHS, IPARM_, &MSGLVL_, &B, &X,
234                              &ERROR);
235    if (ERROR==-4) {
236      // I think this means that the matrix is singular
237      return SYMSOLVER_SINGULAR;
238    }
239    else if (ERROR!=0 ) {
240      Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
241                     "Error in Pardiso during factorization phase.  ERROR = %d.\n", ERROR);
242      return SYMSOLVER_FATAL_ERROR;
243    }
244
245    /* ToDo ask Olaf what this means
246    if (IPARM_[13] != 0) {
247      Jnlst().Printf(J_WARNING, J_LINEAR_ALGEBRA,
248                     "Number of perturbed pivots in factorization phase = %d.\n", IPARM_[13]);
249    }
250    */
251
252    negevals_ = IPARM_[22];
253
254    DBG_ASSERT(IPARM_[21]+IPARM_[22] == dim_);
255
256    // Check whether the number of negative eigenvalues matches the requested
257    // count
258    if (check_NegEVals && (numberOfNegEVals!=negevals_)) {
259      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
260                     "Wrong inertia: required are %d, but we got %d.\n",
261                     numberOfNegEVals, negevals_);
262      return SYMSOLVER_WRONG_INERTIA;
263    }
264
265    return SYMSOLVER_SUCCESS;
266  }
267
268  ESymSolverStatus PardisoSolverInterface::Solve(const Index* ia,
269      const Index* ja,
270      Index nrhs,
271      double *rhs_vals)
272  {
273    DBG_START_METH("PardisoSolverInterface::Solve",dbg_verbosity);
274
275    // Call Pardiso to do the solve for the given right-hand sides
276    ipfint PHASE = 33;
277    ipfint N = dim_;
278    ipfint PERM;   // This should not be accessed by Pardiso
279    ipfint NRHS = nrhs;
280    double* X = new double[nrhs*dim_];
281    ipfint ERROR;
282
283    F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_,
284                              &PHASE, &N, a_, ia, ja, &PERM,
285                              &NRHS, IPARM_, &MSGLVL_, rhs_vals, X,
286                              &ERROR);
287    if (ERROR!=0 ) {
288      Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
289                     "Error in Pardiso during solve phase.  ERROR = %d.\n", ERROR);
290      return SYMSOLVER_FATAL_ERROR;
291    }
292    return SYMSOLVER_SUCCESS;
293  }
294
295  Index PardisoSolverInterface::NumberOfNegEVals() const
296  {
297    DBG_START_METH("PardisoSolverInterface::NumberOfNegEVals",dbg_verbosity);
298    DBG_ASSERT(negevals_>=0);
299    return negevals_;
300  }
301
302  bool PardisoSolverInterface::IncreaseQuality()
303  {
304    // At the moment, I don't see how we could tell Pardiso to do better
305    // (maybe switch from IPARM[20]=1 to IPARM[20]=2?)
306    return false;
307  }
308
309} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.