source: branches/dev/Algorithm/LinearSolvers/IpMa27TSolverInterface.cpp @ 529

Last change on this file since 529 was 529, checked in by andreasw, 15 years ago
  • in Vector class, copy cached values in Copy and update them in Scal
  • perform iterative refinement only once for adaptive strategy
  • fix bugs in PDPerturbationHandler
  • avoid some overhead in CalculateQualityFunction?
  • minor changes in timing
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.8 KB
Line 
1// Copyright (C) 2004, 2005 International Business Machines and others.
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// $Id: IpMa27TSolverInterface.cpp 529 2005-09-29 21:12:38Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2005-03-17
8
9#include "IpMa27TSolverInterface.hpp"
10
11#ifdef HAVE_CMATH
12# include <cmath>
13#else
14# ifdef HAVE_MATH_H
15#  include <math.h>
16# else
17#  error "don't have header file for math"
18# endif
19#endif
20
21/** Prototypes for MA27's Fortran subroutines */
22extern "C"
23{
24  void F77_FUNC(ma27id,MA27ID)(ipfint* ICNTL, double* CNTL);
25  void F77_FUNC(ma27ad,MA27AD)(ipfint *N, ipfint *NZ, const ipfint *IRN, const ipfint* ICN,
26                               ipfint *IW, ipfint* LIW, ipfint* IKEEP, ipfint *IW1,
27                               ipfint* NSTEPS, ipfint* IFLAG, ipfint* ICNTL,
28                               double* CNTL, ipfint *INFO, double* OPS);
29  void F77_FUNC(ma27bd,MA27BD)(ipfint *N, ipfint *NZ, const ipfint *IRN, const ipfint* ICN,
30                               double* A, ipfint* LA, ipfint* IW, ipfint* LIW,
31                               ipfint* IKEEP, ipfint* NSTEPS, ipfint* MAXFRT,
32                               ipfint* IW1, ipfint* ICNTL, double* CNTL,
33                               ipfint* INFO);
34  void F77_FUNC(ma27cd,MA27CD)(ipfint *N, double* A, ipfint* LA, ipfint* IW,
35                               ipfint* LIW, double* W, ipfint* MAXFRT,
36                               double* RHS, ipfint* IW1, ipfint* NSTEPS,
37                               ipfint* ICNTL, double* CNTL);
38}
39
40namespace Ipopt
41{
42#ifdef IP_DEBUG
43  static const Index dbg_verbosity = 0;
44#endif
45
46  Ma27TSolverInterface::Ma27TSolverInterface()
47      :
48      dim_(0),
49      nonzeros_(0),
50      initialized_(false),
51      pivtol_changed_(false),
52      refactorize_(false),
53
54      liw_(0),
55      iw_(NULL),
56      ikeep_(NULL),
57      la_(0),
58      a_(NULL),
59
60      la_increase_(false),
61      liw_increase_(false)
62  {
63    DBG_START_METH("Ma27TSolverInterface::Ma27TSolverInterface()",dbg_verbosity);
64  }
65
66  Ma27TSolverInterface::~Ma27TSolverInterface()
67  {
68    DBG_START_METH("Ma27TSolverInterface::~Ma27TSolverInterface()",
69                   dbg_verbosity);
70    delete [] iw_;
71    delete [] ikeep_;
72    delete [] a_;
73  }
74
75  void Ma27TSolverInterface::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
76  {
77    roptions->AddBoundedNumberOption(
78      "pivtol",
79      "Pivot tolerance for the linear solver MA27.",
80      0.0, true, 1.0, true, 1e-8,
81      "A smaller number pivots for sparsity, "
82      "a larger number pivots for stability.");
83    roptions->AddBoundedNumberOption(
84      "pivtolmax",
85      "Maximum pivot tolerance.",
86      0.0, true, 1.0, true, 1e-4,
87      "Ipopt may increase pivtol as high as pivtolmax "
88      "to get a more accurate solution to the linear system.");
89    roptions->AddLowerBoundedNumberOption(
90      "liw_init_factor",
91      "Integer workspace memory for MA27.",
92      1.0, false, 5.0,
93      "The initial integer workspace memory = liw_init_factor * memory "
94      "required by unfactored system. Ipopt will increase the workspace "
95      "size by meminc_factor if required.");
96    roptions->AddLowerBoundedNumberOption(
97      "la_init_factor",
98      "Real workspace memory for MA27.",
99      1.0, false, 5.0,
100      "The initial real workspace memory = la_init_factor * memory "
101      "required by unfactored system. Ipopt will increase the workspace"
102      " size by meminc_factor if required.");
103    roptions->AddLowerBoundedNumberOption(
104      "meminc_factor",
105      "Increment factor for workspace size for MA27.",
106      1.0, false, 10.0,
107      "If the integer or real workspace is not large enough, "
108      "Ipopt will increase its size by this factor.");
109  }
110
111  bool Ma27TSolverInterface::InitializeImpl(const OptionsList& options,
112      const std::string& prefix)
113  {
114    options.GetNumericValue("pivtol", pivtol_, prefix);
115    if(options.GetNumericValue("pivtolmax", pivtolmax_, prefix)) {
116      ASSERT_EXCEPTION(pivtolmax_>=pivtol_, OPTION_INVALID,
117                       "Option \"pivtolmax\": This value must be between "
118                       "pivtol and 1.");
119    }
120    else {
121      pivtolmax_ = Max(pivtolmax_, pivtol_);
122    }
123
124    options.GetNumericValue("liw_init_factor", liw_init_factor_, prefix);
125    options.GetNumericValue("la_init_factor", la_init_factor_, prefix);
126    options.GetNumericValue("meminc_factor", meminc_factor_, prefix);
127    // The following option is registered by OrigIpoptNLP
128    options.GetBoolValue("warm_start_same_structure",
129                         warm_start_same_structure_, prefix);
130
131    /* Set the default options for MA27 */
132    F77_FUNC(ma27id,MA27ID)(icntl_, cntl_);
133#ifndef IP_DEBUG
134
135    icntl_[0] = 0;       // Suppress error messages
136    icntl_[1] = 0;       // Suppress diagnostic messages
137#endif
138
139    // Reset all private data
140    initialized_=false;
141    pivtol_changed_ = false;
142    refactorize_ = false;
143
144    la_increase_=false;
145    liw_increase_=false;
146
147    if (!warm_start_same_structure_) {
148      dim_=0;
149      nonzeros_=0;
150    }
151    else {
152      ASSERT_EXCEPTION(dim_>0 && nonzeros_>0, INVALID_WARMSTART,
153                       "Ma27TSolverInterface called with warm_start_same_structure, but the problem is solved for the first time.");
154    }
155
156    return true;
157  }
158
159  ESymSolverStatus Ma27TSolverInterface::MultiSolve(bool new_matrix,
160      const Index* airn,
161      const Index* ajcn,
162      Index nrhs,
163      double* rhs_vals,
164      bool check_NegEVals,
165      Index numberOfNegEVals)
166  {
167    DBG_START_METH("Ma27TSolverInterface::MultiSolve",dbg_verbosity);
168    DBG_ASSERT(!check_NegEVals || ProvidesInertia());
169    DBG_ASSERT(initialized_);
170    DBG_ASSERT(la_!=0);
171
172    if (pivtol_changed_) {
173      DBG_PRINT((1,"Pivot tolerance has changed.\n"));
174      pivtol_changed_ = false;
175      // If the pivot tolerance has been changed but the matrix is not
176      // new, we have to request the values for the matrix again to do
177      // the factorization again.
178      if (!new_matrix) {
179        DBG_PRINT((1,"Ask caller to call again.\n"));
180        refactorize_ = true;
181        return SYMSOLVER_CALL_AGAIN;
182      }
183    }
184
185    // check if a factorization has to be done
186    DBG_PRINT((1, "new_matrix = %d\n", new_matrix));
187    if (new_matrix || refactorize_) {
188      // perform the factorization
189      ESymSolverStatus retval;
190      retval = Factorization(airn, ajcn, check_NegEVals, numberOfNegEVals);
191      if (retval!=SYMSOLVER_SUCCESS) {
192        DBG_PRINT((1, "FACTORIZATION FAILED!\n"));
193        return retval;  // Matrix singular or error occurred
194      }
195      refactorize_ = false;
196    }
197
198    // do the backsolve
199    return Backsolve(nrhs, rhs_vals);
200  }
201
202  double* Ma27TSolverInterface::GetValuesArrayPtr()
203  {
204    DBG_START_METH("Ma27TSolverInterface::GetValuesArrayPtr",dbg_verbosity);
205    DBG_ASSERT(initialized_);
206
207    // If the size of a is to be increase for the next factorization
208    // anyway, delete the current large array and just return enough
209    // to store the values
210
211    if (la_increase_) {
212      delete [] a_;
213      a_ = new double [nonzeros_];
214    }
215
216    return a_;
217  }
218
219  /** Initialize the local copy of the positions of the nonzero
220      elements */
221  ESymSolverStatus Ma27TSolverInterface::InitializeStructure(Index dim, Index nonzeros,
222      const Index* airn,
223      const Index* ajcn)
224  {
225    DBG_START_METH("Ma27TSolverInterface::InitializeStructure",dbg_verbosity);
226
227    ESymSolverStatus retval = SYMSOLVER_SUCCESS;
228    if (!warm_start_same_structure_) {
229      dim_ = dim;
230      nonzeros_ = nonzeros;
231
232      // Do the symbolic facotrization
233      retval = SymbolicFactorization(airn, ajcn);
234      if (retval != SYMSOLVER_SUCCESS ) {
235        return retval;
236      }
237    }
238    else {
239      ASSERT_EXCEPTION(dim_==dim && nonzeros_==nonzeros, INVALID_WARMSTART,
240                       "Ma27TSolverInterface called with warm_start_same_structure, but the problem size has changed.");
241    }
242
243    initialized_ = true;
244
245    return retval;
246  }
247
248  ESymSolverStatus Ma27TSolverInterface::SymbolicFactorization(const Index* airn,
249      const Index* ajcn)
250  {
251    DBG_START_METH("Ma27TSolverInterface::SymbolicFactorization",dbg_verbosity);
252
253    IpData().TimingStats().LinearSystemSymbolicFactorization().Start();
254
255    // Get memory for the IW workspace
256    delete [] iw_;
257
258    // Overstimation factor for LIW (20% recommended in MA27 documentation)
259    const double LiwFact = 2.0;   // This is 100% overestimation
260    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
261                   "In Ma27TSolverInterface::InitializeStructure: Using overestimation factor LiwFact = %e\n",
262                   LiwFact);
263    liw_ = (ipfint)(LiwFact*(double(2*nonzeros_+3*dim_+1)));
264    iw_ = new ipfint[liw_];
265
266    // Get memory for IKEEP
267    delete [] ikeep_;
268    ikeep_ = new ipfint[3*dim_];
269
270    // Call MA27AD (cast to ipfint for Index types)
271    ipfint N = dim_;
272    ipfint NZ = nonzeros_;
273    ipfint IFLAG = 0;
274    double OPS;
275    ipfint INFO[20];
276    ipfint* IW1 = new ipfint[2*dim_];  // Get memory for IW1 (only local)
277    F77_FUNC(ma27ad,MA27AD)(&N, &NZ, airn, ajcn, iw_, &liw_, ikeep_,
278                            IW1, &nsteps_, &IFLAG, icntl_, cntl_,
279                            INFO, &OPS);
280    delete [] IW1;  // No longer required
281
282    // Receive several information
283    ipfint iflag = INFO[0];   // Information flag
284    ipfint ierror = INFO[1];  // Error flag
285    ipfint nrlnec = INFO[4];  // recommended value for la
286    ipfint nirnec = INFO[5];  // recommended value for liw
287
288    // Check if error occurred
289    if (iflag!=0) {
290      Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
291                     "*** Error from MA27AD *** IFLAG = %d IERROR = %d\n", iflag, ierror);
292      if (iflag==1) {
293        Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
294                       "The index a matrix is out of range.\nPlease check your implementation of the Jabobian and Hessian matrices.");
295      }
296      return SYMSOLVER_FATAL_ERROR;
297    }
298
299    // ToDo: try and catch
300    // Reserve memory for iw_ for later calls, based on suggested size
301    delete [] iw_;
302    liw_ = (ipfint)(liw_init_factor_ * (double)(nirnec));
303    iw_ = new ipfint[liw_];
304
305    // Reserve memory for a_
306    delete [] a_;
307    la_ = Max(nonzeros_,(ipfint)(la_init_factor_ * (double)(nrlnec)));
308    a_ = new double[la_];
309
310    IpData().TimingStats().LinearSystemSymbolicFactorization().End();
311
312    return SYMSOLVER_SUCCESS;
313  }
314
315  ESymSolverStatus
316  Ma27TSolverInterface::Factorization(const Index* airn,
317                                      const Index* ajcn,
318                                      bool check_NegEVals,
319                                      Index numberOfNegEVals)
320  {
321    DBG_START_METH("Ma27TSolverInterface::Factorization",dbg_verbosity);
322    // Check if la should be increased
323    IpData().TimingStats().LinearSystemFactorization().Start();
324    if (la_increase_) {
325      double* a_old = a_;
326      ipfint la_old = la_;
327      la_ = (ipfint)(meminc_factor_ * (double)(la_));
328      a_ = new double[la_];
329      for (Index i=0; i<nonzeros_; i++) {
330        a_[i] = a_old[i];
331      }
332      delete [] a_old;
333      la_increase_ = false;
334      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
335                     "In Ma27TSolverInterface::Factorization: Increasing la from %d to %d\n",
336                     la_old, la_);
337    }
338
339    // Check if liw should be increased
340    if (liw_increase_) {
341      delete [] iw_;
342      ipfint liw_old = liw_;
343      liw_ = (ipfint)(meminc_factor_ * (double)(liw_));
344      iw_ = new ipfint[liw_];
345      liw_increase_ = false;
346      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
347                     "In Ma27TSolverInterface::Factorization: Increasing liw from %d to %d\n",
348                     liw_old, liw_);
349    }
350
351    ipfint iflag;  // Information flag
352    ipfint ncmpbr; // Number of double precision compressions
353    ipfint ncmpbi; // Number of integer compressions
354
355    // Call MA27BD; possibly repeatedly if workspaces are too small
356    ipfint N=dim_;
357    ipfint NZ=nonzeros_;
358    ipfint* IW1 = new ipfint[2*dim_];
359    ipfint INFO[20];
360    cntl_[0] = pivtol_;  // Set pivot tolerance
361
362    F77_FUNC(ma27bd,MA27BD)(&N, &NZ, airn, ajcn, a_,
363                            &la_, iw_, &liw_, ikeep_, &nsteps_,
364                            &maxfrt_, IW1, icntl_, cntl_, INFO);
365    delete [] IW1;
366
367    // Receive information about the factorization
368    iflag = INFO[0];        // Information flag
369    ipfint ierror = INFO[1];  // Error flag
370    ncmpbr = INFO[11];      // Number of double compressions
371    ncmpbi = INFO[12];      // Number of integer compressions
372    negevals_ = INFO[14];   // Number of negative eigenvalues
373
374    DBG_PRINT((1,"Return from MA27 iflag = %d and ierror = %d\n",
375               iflag, ierror));
376
377    // Check if factorization failed due to insufficient memory space
378    // iflag==-3 if LIW too small (recommended value in ierror)
379    // iflag==-4 if LA too small (recommended value in ierror)
380    if (iflag==-3 || iflag==-4) {
381      // Increase size of both LIW and LA
382      delete [] iw_;
383      delete [] a_;
384      ipfint liw_old = liw_;
385      ipfint la_old = la_;
386      if(iflag==-3) {
387        liw_ = (ipfint)(meminc_factor_ * (double)(ierror));
388        la_ = (ipfint)(meminc_factor_ * (double)(la_));
389      }
390      else {
391        liw_ = (ipfint)(meminc_factor_ * (double)(liw_));
392        la_ = (ipfint)(meminc_factor_ * (double)(ierror));
393      }
394      iw_ = new ipfint[liw_];
395      a_ = new double[la_];
396      Jnlst().Printf(J_WARNING, J_LINEAR_ALGEBRA,
397                     "MA27BD returned iflag=%d.\n Increase liw from %d to %d and la from %d to %d and factorize again.\n",
398                     iflag, liw_old, liw_, la_old, la_);
399      IpData().TimingStats().LinearSystemFactorization().End();
400      return SYMSOLVER_CALL_AGAIN;
401    }
402
403    // Check if the system is singular, and if some other error occurred
404    if (iflag==-5 || iflag==3) {
405      IpData().TimingStats().LinearSystemFactorization().End();
406      return SYMSOLVER_SINGULAR;
407    }
408    else if (iflag != 0) {
409      // There is some error
410      IpData().TimingStats().LinearSystemFactorization().End();
411      return SYMSOLVER_FATAL_ERROR;
412    }
413
414    // Check if it might be more efficient to use more memory next time
415    // (if there were too many compressions for this factorization)
416    if (ncmpbr>=10) {
417      la_increase_ = true;
418      Jnlst().Printf(J_WARNING, J_LINEAR_ALGEBRA,
419                     "MA27BD returned ncmpbr=%d. Increase la before the next factorization.\n",
420                     ncmpbr);
421    }
422    if (ncmpbi>=10) {
423      liw_increase_ = true;
424      Jnlst().Printf(J_WARNING, J_LINEAR_ALGEBRA,
425                     "MA27BD returned ncmpbi=%d. Increase liw before the next factorization.\n",
426                     ncmpbr);
427    }
428
429    // Check whether the number of negative eigenvalues matches the requested
430    // count
431    if (check_NegEVals && (numberOfNegEVals!=negevals_)) {
432      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
433                     "In Ma27TSolverInterface::Factorization: negevals_ = %d, but numberOfNegEVals = %d\n",
434                     negevals_, numberOfNegEVals);
435      IpData().TimingStats().LinearSystemFactorization().End();
436      return SYMSOLVER_WRONG_INERTIA;
437    }
438
439    IpData().TimingStats().LinearSystemFactorization().End();
440    return SYMSOLVER_SUCCESS;
441  }
442
443  ESymSolverStatus Ma27TSolverInterface::Backsolve(Index nrhs,
444      double *rhs_vals)
445  {
446    DBG_START_METH("Ma27TSolverInterface::Backsolve",dbg_verbosity);
447    IpData().TimingStats().LinearSystemBackSolve().Start();
448
449    ipfint N=dim_;
450    double* W = new double[maxfrt_];
451    ipfint* IW1 = new ipfint[nsteps_];
452
453    // For each right hand side, call MA27CD
454    for(Index irhs=0; irhs<nrhs; irhs++) {
455      if (DBG_VERBOSITY()>=2) {
456        for (Index i=0; i<dim_; i++) {
457          DBG_PRINT((2, "rhs[%5d] = %23.15e\n", i, rhs_vals[irhs*dim_+i]));
458        }
459      }
460
461      F77_FUNC(ma27cd,MA27CD)(&N, a_, &la_, iw_, &liw_, W, &maxfrt_,
462                              &rhs_vals[irhs*dim_], IW1, &nsteps_,
463                              icntl_, cntl_);
464
465      if (DBG_VERBOSITY()>=2) {
466        for (Index i=0; i<dim_; i++) {
467          DBG_PRINT((2, "sol[%5d] = %23.15e\n", i, rhs_vals[irhs*dim_+i]));
468        }
469      }
470    }
471    delete [] W;
472    delete [] IW1;
473
474    IpData().TimingStats().LinearSystemBackSolve().End();
475    return SYMSOLVER_SUCCESS;
476  }
477
478  Index Ma27TSolverInterface::NumberOfNegEVals() const
479  {
480    DBG_START_METH("Ma27TSolverInterface::NumberOfNegEVals",dbg_verbosity);
481    DBG_ASSERT(ProvidesInertia());
482    DBG_ASSERT(initialized_);
483    return negevals_;
484  }
485
486  bool Ma27TSolverInterface::IncreaseQuality()
487  {
488    DBG_START_METH("Ma27TSolverInterface::IncreaseQuality",dbg_verbosity);
489    if (pivtol_ == pivtolmax_) {
490      return false;
491    }
492    pivtol_changed_ = true;
493
494    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
495                   "Indreasing pivot tolerance from %7.2e ",
496                   pivtol_);
497    pivtol_ = Min(pivtolmax_, pow(pivtol_,0.75));
498    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
499                   "to %7.2e.\n",
500                   pivtol_);
501    return true;
502  }
503
504} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.