source: branches/dev/Algorithm/LinearSolvers/IpMa57TSolverInterface.cpp @ 566

Last change on this file since 566 was 566, checked in by andreasw, 15 years ago

astyle

  • Property svn:eol-style set to native
  • Property svn:keywords set to "Author Date Id Revision"
File size: 24.9 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: IpMa57TSolverInterface.cpp 566 2005-11-01 18:09:34Z andreasw $
6//
7// Authors:  Michael Hagemann               Univ of Basel 2005-10-28
8//               original version (based on MA27TSolverInterface.cpp)
9
10#include "IpMa57TSolverInterface.hpp"
11
12#ifdef HAVE_CMATH
13# include <cmath>
14#else
15# ifdef HAVE_MATH_H
16#  include <math.h>
17# else
18#  error "don't have header file for math"
19# endif
20#endif
21
22#include <iostream>
23
24
25/** Prototypes for MA57's Fortran subroutines */
26extern "C"
27{
28  /*
29   *  MA57ID -- Initialize solver.
30   */
31  extern void  F77_FUNC (ma57id, MA57ID) (
32      double    *cntl,
33      ipfint    *icntl);
34
35  /*
36   *  MA57AD -- Symbolic Factorization.
37   */
38  extern void  F77_FUNC (ma57ad, MA57AD) (
39      ipfint    *n,     /* Order of matrix. */
40      ipfint    *ne,            /* Number of entries. */
41
42      const ipfint    *irn,       /* Matrix nonzero row structure */
43      const ipfint    *jcn,       /* Matrix nonzero column structure */
44
45      ipfint    *lkeep,     /* Workspace for the pivot order of lenght 3*n */
46      ipfint    *keep,      /* Workspace for the pivot order of lenght 3*n */
47      /* Automatically iflag = 0; ikeep pivot order iflag = 1 */
48      ipfint    *iwork,     /* Integer work space. */
49      ipfint    *icntl,     /* Integer Control parameter of length 30*/
50      ipfint    *info,      /* Statistical Information; Integer array of length 20 */
51      double    *rinfo);    /* Double Control parameter of length 5 */
52
53  /*
54   * MA57BD -- Numerical Factorization.
55   */
56  extern void  F77_FUNC (ma57bd, MA57BD) (
57      ipfint    *n,     /* Order of matrix. */
58      ipfint    *ne,            /* Number of entries. */
59
60      double    *a,     /* Numerical values. */
61      double    *fact,      /* Entries of factors. */
62      ipfint    *lfact,     /* Length of array `fact'. */
63      ipfint    *ifact,     /* Indexing info for factors. */
64      ipfint    *lifact,    /* Length of array `ifact'. */
65
66      ipfint    *lkeep,     /* Length of array `keep'. */
67      ipfint    *keep,      /* Integer array. */
68
69      ipfint    *iwork,     /* Workspace of length `n'. */
70
71      ipfint    *icntl,     /* Integer Control parameter of length 20. */
72      double    *cntl,      /* Double Control parameter of length 5. */
73      ipfint    *info,      /* Statistical Information; Integer array of length 40. */
74      double    *rinfo);    /* Statistical Information; Real array of length 20. */
75
76  /*
77   * MA57CD -- Solution.
78   */
79  extern void  F77_FUNC (ma57cd, MA57CD) (
80      ipfint    *job,       /* Solution job.  Solve for...
81                                                                   JOB <= 1:  A
82                                                                   JOB == 2:  PLP^t
83                                                                   JOB == 3:  PDP^t
84                                                                   JOB >= 4:  PL^t P^t */
85
86      ipfint    *n,     /* Order of matrix. */
87
88      double    *fact,      /* Entries of factors. */
89      ipfint    *lfact,     /* Length of array `fact'. */
90      ipfint    *ifact,     /* Indexing info for factors. */
91      ipfint    *lifact,    /* Length of array `ifact'. */
92
93      ipfint    *nrhs,      /* Number of right hand sides. */
94      double    *rhs,       /* Numerical Values. */
95      ipfint    *lrhs,      /* Leading dimensions of `rhs'. */
96
97      double    *work,      /* Real workspace. */
98      ipfint    *lwork,     /* Length of `work', >= N*NRHS. */
99      ipfint    *iwork,     /* Integer array of length `n'. */
100
101      ipfint    *icntl,     /* Integer Control parameter array of length 20. */
102      ipfint    *info);     /* Statistical Information; Integer array of length 40. */
103
104  /*
105   * MC57ED -- Copy arrays.
106   */
107  extern void  F77_FUNC (ma57ed, MA57ED) (
108      ipfint    *n,
109      ipfint    *ic,        /* 0: copy real array.  >=1:  copy integer array. */
110      ipfint    *keep,
111
112      double    *fact,
113      ipfint    *lfact,
114      double    *newfac,
115      ipfint    *lnew,
116
117      ipfint    *ifact,
118      ipfint    *lifact,
119      ipfint    *newifc,
120      ipfint    *linew,
121
122      ipfint    *info);
123}
124
125namespace Ipopt
126{
127#ifdef IP_DEBUG
128  static const Index dbg_verbosity = 0;
129#endif
130
131  std::string ma57_err_msg[] = {
132                                 "Operation successful.\n",
133
134                                 "Value of N is out of range on a call to MA57A/AD, MA57B/BD, MA57C/CD, or\n"
135                                 "MA57D/DD. Value given is held in INFO(2).\n",
136
137                                 "Value of NE is out of range on a call to MA57A/AD, MA57B/BD, or\n"
138                                 "MA57D/DD. Value given is held in INFO(2).\n",
139
140                                 "Failure due to insufficient REAL space on a call to MA57B/BD. INFO(17)\n"
141                                 "is set to a value that may suffice. INFO(2) is set to value of\n"
142                                 "LFACT. The user can allocate a larger array and copy the contents of\n"
143                                 "FACT into it using MA57E/ED, and recall MA57B/BD.\n",
144
145                                 "Failure due to insufficient INTEGER space on a call to\n"
146                                 "MA57B/BD. INFO(18) is set to a value that may suffice. INFO(2) is set to\n"
147                                 "value of LIFACT. The user can allocate a larger array and copy the\n"
148                                 "contents of IFACT into it using MA57E/ED, and recall MA57B/BD.\n",
149
150                                 "A pivot with magnitude less than or equal to CNTL(2) was found at pivot\n"
151                                 "step INFO(2) when calling MA57B/BD with ICNTL(7) = 2 or 3, or the\n"
152                                 "correction obtained when using matrix modification does not give a pivot\n"
153                                 "greater than CNTL(2) when ICNTL(7) = 4.\n",
154
155                                 "A change of sign of pivots has been detected when ICNTL(7) = 2. INFO(2)\n"
156                                 "is set to the pivot step at which the change was detected on a call to\n"
157                                 "MA57B/BD.\n",
158
159                                 "Either LNEW < LFACT or LINEW < LIFACT on a call to MA57E/ED. INFO(2) is\n"
160                                 "set to LNEW or LINEW as appropriate.\n",
161
162                                 "Iterative refinement fails to converge in specified number of iterations\n"
163                                 "on a call to MA57D/DD.\n",
164
165                                 "Error in permutation array when ICNTL(6)=1 on a call to\n"
166                                 "MA57A/AD. INFO(2) holds first component at which error was detected.\n",
167
168                                 "Value of ICNTL(7) out of range on a call to MA57B/BD. Value given held\n"
169                                 "in INFO(2).\n",
170
171                                 "LRHS < N on a call to MA57C/CD. INFO(2) holds value of LRHS.\n",
172
173                                 "Invalid value for JOB on a call to MA57D/DD. Value given held in\n"
174                                 "INFO(2).\n",
175
176                                 "Invalid value of ICNTL(9) on a call to MA57D/DD. Value given held in\n"
177                                 "INFO(2).\n",
178
179                                 "Failure of MC71A/AD on a call to MA57D/DD with ICNTL(10)> 0.\n",
180
181                                 "LKEEP less than 5*N+NE+MAX(N,NE) +42 on a call to MA57A/AD or\n"
182                                 "MA57B/BD. INFO(2) holds value of LKEEP.\n",
183
184                                 "NRHS less than 1 on call to MA57C/CD. INFO(2) holds value of NRHS.\n",
185
186                                 "LWORK too small on entry to MA57C/CD. INFO(2) holds minimum length\n"
187                                 "required. A positive value of INFO(1) is associated with a warning\n"
188                                 "message that will be output on unit ICNTL(2).\n"
189                               };
190
191  std::string ma57_wrn_msg[] = {
192                                 "Operation successful.\n",
193
194                                 "Index (in IRN or JCN) out of range on call to MA57A/AD or\n"
195                                 "MA57D/DD. Action taken by subroutine is to ignore any such entries and\n"
196                                 "continue. INFO(3) is set to the number of faulty entries. Details of the\n"
197                                 "first ten are printed on unit ICNTL(2).\n",
198
199                                 "Duplicate indices on call to MA57A/AD or MA57D/DD. Action taken by\n"
200                                 "subroutine is to keep the duplicates and then to sum corresponding reals\n"
201                                 "when MA57B/BD is called. INFO(4) is set to the number of faulty\n"
202                                 "entries. Details of the first ten are printed on unit ICNTL(2).\n",
203
204                                 "Both out-of-range indices and duplicates exist.\n",
205
206                                 "Matrix is rank deficient on exit from MA57B/BD. In this case, a\n"
207                                 "decomposition will still have been produced that will enable the\n"
208                                 "subsequent solution of consistent equations. INFO(25) will be set to the\n"
209                                 "rank of the factorized matrix.\n",
210
211                                 "Pivots have different signs when factorizing a supposedly definite\n"
212                                 "matrix (ICNTL(7) = 3) on call to MA57B/BD. INFO(26) is set to the number\n"
213                                 "of sign changes.\n",
214
215                                 "-",
216                                 "-",
217
218                                 "During error analysis the infinity norm of the computed solution was\n"
219                                 "found to be zero.\n",
220
221                                 "Insufficient real space to complete factorization when MA57B/BD called\n"
222                                 "with ICNTL(8) != 0. User can copy real values to a longer array using\n"
223                                 "MA57E/ED and recall MA57B/BD using this longer array to continue the\n"
224                                 "factorization.\n",
225
226                                 "Insufficient integer space to complete factorization when MA57B/BD\n"
227                                 "called with ICNTL(8) != 0. User can copy integer values to a longer\n"
228                                 "array using MA57E/ED and recall MA57B/BD using this longer array to\n"
229                                 "continue the factorization.\n"
230                               };
231
232  Ma57TSolverInterface::Ma57TSolverInterface()
233      :
234      dim_(0),
235      nonzeros_(0),
236      initialized_(false),
237      pivtol_changed_(false),
238      refactorize_(false),
239
240      pivtol_(0.01),
241      warm_start_same_structure_(false),
242      wd_keep_(NULL),
243      wd_iwork_(NULL),
244      wd_fact_(NULL),
245      wd_ifact_(NULL),
246      a_(NULL)
247  {
248    DBG_START_METH("Ma57TSolverInterface::Ma57TSolverInterface()",
249                   dbg_verbosity);
250  }
251
252  Ma57TSolverInterface::~Ma57TSolverInterface()
253  {
254    DBG_START_METH("Ma57TSolverInterface::~Ma57TSolverInterface()",
255                   dbg_verbosity);
256    delete [] a_;
257
258    delete [] wd_fact_;
259    delete [] wd_ifact_;
260
261    delete [] wd_iwork_;
262    delete [] wd_keep_;
263  }
264
265  void Ma57TSolverInterface::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
266  {
267    roptions->AddBoundedNumberOption(
268      "ma57_pivtol",
269      "Pivot tolerance for the linear solver MA57.",
270      0.0, true, 1.0, true, 1e-8,
271      "A smaller number pivots for sparsity, "
272      "a larger number pivots for stability.");
273    roptions->AddBoundedNumberOption(
274      "ma57_pivtolmax",
275      "Maximum pivot tolerance for the linear solver MA57.",
276      0.0, true, 1.0, true, 1e-4,
277      "Ipopt may increase pivtol as high as pivtolmax "
278      "to get a more accurate solution to the linear system.");
279    roptions->AddLowerBoundedNumberOption(
280      "ma57_pre_alloc",
281      "Safety factor for work space memory allocation for the linear solver MA57.",
282      1., false, 3.,
283      "If 1 is chosen, the suggested amount of work space is used.  However, "
284      "choosen a larger number might avoid reallocation if the suggest values "
285      "do not suffice.");
286  }
287
288  bool Ma57TSolverInterface::InitializeImpl(const OptionsList&  options,
289      const std::string&    prefix)
290  {
291    // Obtain the options settings
292    options.GetNumericValue("ma57_pivtol", pivtol_, prefix);
293    if(options.GetNumericValue("ma57_pivtolmax", pivtolmax_, prefix)) {
294      ASSERT_EXCEPTION(pivtolmax_>=pivtol_, OPTION_INVALID,
295                       "Option \"pivtolmax\": This value must be between "
296                       "pivtol and 1.");
297    }
298    else {
299      pivtolmax_ = Max(pivtolmax_, pivtol_);
300    }
301
302    options.GetNumericValue("ma57_pre_alloc", ma57_pre_alloc_, prefix);
303
304    // The following option is registered by OrigIpoptNLP
305    options.GetBoolValue("warm_start_same_structure",
306                         warm_start_same_structure_, prefix);
307    DBG_ASSERT(!warm_start_same_structure_ && "warm_start_same_structure not yet implemented");
308
309
310    /* Initialize. */
311    F77_FUNC (ma57id, MA57ID) (wd_cntl_, wd_icntl_);
312
313    /* Custom settings for MA57. */
314    wd_icntl_[1-1] = 0;      /* Error stream */
315    wd_icntl_[2-1] = 0;      /* Warning stream. */
316
317    wd_icntl_[4-1] = 1;      /* Print statistics.  NOT Used. */
318    wd_icntl_[5-1] = 0;      /* Print error. */
319
320    // wd_icntl[6-1] = 0;       /* Pivoting order. */
321
322    wd_cntl_[1-1]  = pivtol_;    /* Pivot threshold. */
323    wd_icntl_[7-1] = 1;      /* Pivoting strategy. */
324
325    // wd_icntl[8-1] = 0;       /* Retry factorization. */
326
327    if (!warm_start_same_structure_) {
328      dim_=0;
329      nonzeros_=0;
330      delete [] a_;
331      delete [] wd_fact_;
332      delete [] wd_ifact_;
333      delete [] wd_iwork_;
334      delete [] wd_keep_;
335    }
336    else {
337      ASSERT_EXCEPTION(dim_>0 && nonzeros_>0, INVALID_WARMSTART,
338                       "Ma57TSolverInterface called with warm_start_same_structure, "
339                       "but the problem is solved for the first time.");
340    }
341
342    return true;
343  }
344
345  ESymSolverStatus
346  Ma57TSolverInterface::MultiSolve(bool         new_matrix,
347                                   const Index*     airn,
348                                   const Index*     ajcn,
349                                   Index        nrhs,
350                                   double*      rhs_vals,
351                                   bool         check_NegEVals,
352                                   Index        numberOfNegEVals)
353  {
354    DBG_START_METH("Ma57TSolverInterface::MultiSolve",dbg_verbosity);
355
356    // DBG_ASSERT(!check_NegEVals || ProvidesInertia());
357    // DBG_ASSERT(initialized_);
358    // DBG_ASSERT(la_!=0);
359
360    if (pivtol_changed_) {
361      DBG_PRINT((1,"Pivot tolerance has changed.\n"));
362      pivtol_changed_ = false;
363      // If the pivot tolerance has been changed but the matrix is not
364      // new, we have to request the values for the matrix again to do
365      // the factorization again.
366      if (!new_matrix) {
367        DBG_PRINT((1,"Ask caller to call again.\n"));
368        refactorize_ = true;
369        return SYMSOLVER_CALL_AGAIN;
370      }
371    }
372
373    // check if a factorization has to be done
374    DBG_PRINT((1, "new_matrix = %d\n", new_matrix));
375    if (new_matrix || refactorize_) {
376      // perform the factorization
377      ESymSolverStatus retval;
378      retval = Factorization(airn, ajcn, check_NegEVals, numberOfNegEVals);
379      if (retval!=SYMSOLVER_SUCCESS) {
380        DBG_PRINT((1, "FACTORIZATION FAILED!\n"));
381        return retval;  // Matrix singular or error occurred
382      }
383      refactorize_ = false;
384    }
385
386    // do the backsolve
387    return Backsolve(nrhs, rhs_vals);
388  }
389
390  double* Ma57TSolverInterface::GetValuesArrayPtr()
391  {
392    DBG_START_METH("Ma57TSolverInterface::GetValuesArrayPtr",dbg_verbosity);
393    DBG_ASSERT(initialized_);
394
395    return a_;
396  }
397
398  /** Initialize the local copy of the positions of the nonzero
399      elements */
400  ESymSolverStatus
401  Ma57TSolverInterface::InitializeStructure(
402    Index       dim,
403    Index     nonzeros,
404    const Index*  airn,
405    const Index*  ajcn)
406  {
407    DBG_START_METH("Ma57TSolverInterface::InitializeStructure",dbg_verbosity);
408
409    ESymSolverStatus retval = SYMSOLVER_SUCCESS;
410    if (!warm_start_same_structure_) {
411      dim_ = dim;
412      nonzeros_ = nonzeros;
413      // for MA57, a_ only has to be as long as the number of nonzero
414      // elements
415      delete [] a_;
416      a_ = new double [nonzeros_];
417
418      // Do the symbolic facotrization
419      retval = SymbolicFactorization(airn, ajcn);
420      if (retval != SYMSOLVER_SUCCESS ) {
421        return retval;
422      }
423    }
424    else {
425      ASSERT_EXCEPTION(dim_==dim && nonzeros_==nonzeros, INVALID_WARMSTART,
426                       "Ma57TSolverInterface called with warm_start_same_structure, "
427                       "but the problem size has changed.");
428    }
429
430    initialized_ = true;
431
432    return retval;
433  }
434
435  ESymSolverStatus
436  Ma57TSolverInterface::SymbolicFactorization(const Index*  airn,
437      const Index*  ajcn)
438  {
439    DBG_START_METH("Ma57TSolverInterface::SymbolicFactorization",dbg_verbosity);
440    IpData().TimingStats().LinearSystemSymbolicFactorization().Start();
441
442    ipfint n  = dim_;
443    ipfint ne = nonzeros_;
444
445    wd_lkeep_ = 5*n + ne + Max (n,ne) + 42;
446
447    wd_iwork_ = new ipfint[5*n];
448    wd_keep_  = new ipfint[wd_lkeep_];
449
450    F77_FUNC (ma57ad, MA57AD)
451    (&n, &ne, airn, ajcn, &wd_lkeep_, wd_keep_, wd_iwork_,
452     wd_icntl_, wd_info_, wd_rinfo_);
453
454    if (wd_info_[0] < 0) {
455      Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
456                     "*** Error from MA57AD *** INFO(0) = %d\n", wd_info_[0]);
457    }
458
459    wd_lfact_  = (ipfint)((Number)wd_info_[8] * ma57_pre_alloc_);
460    wd_lifact_ = (ipfint)((Number)wd_info_[9] * ma57_pre_alloc_);
461
462    // XXX MH:  Why is this necessary?  Is `::Factorization' called more
463    // than once per object lifetime?  Where should allocation take
464    // place, then?
465
466    // AW: I moved this here now - my understanding is that wd_info[8]
467    // and wd_info[9] are computed here, so we can just allocate the
468    // amount of memory here.  I don't think there is any need to
469    // reallocate it later for every factorization
470    delete [] wd_fact_;
471    delete [] wd_ifact_;
472
473    wd_fact_  = new double[wd_lfact_];
474    wd_ifact_ = new int[wd_lifact_];
475
476    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
477                   "Suggested lfact  (*%e):  %d\n", ma57_pre_alloc_, wd_lfact_);
478    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
479                   "Suggested lifact (*%e):  %d\n", ma57_pre_alloc_, wd_lifact_);
480
481    IpData().TimingStats().LinearSystemSymbolicFactorization().End();
482    return SYMSOLVER_SUCCESS;
483  }
484
485  ESymSolverStatus
486  Ma57TSolverInterface::Factorization(const Index*  airn,
487                                      const Index*  ajcn,
488                                      bool          check_NegEVals,
489                                      Index         numberOfNegEVals)
490  {
491    DBG_START_METH("Ma57TSolverInterface::Factorization",dbg_verbosity);
492    IpData().TimingStats().LinearSystemFactorization().Start();
493
494    int fact_error = 1;
495
496    ipfint n  = dim_;
497    ipfint ne = nonzeros_;
498
499    while (fact_error > 0) {
500      F77_FUNC (ma57bd, MA57BD)
501      (&n, &ne, a_, wd_fact_, &wd_lfact_, wd_ifact_, &wd_lifact_,
502       &wd_lkeep_, wd_keep_, wd_iwork_,
503       wd_icntl_, wd_cntl_, wd_info_, wd_rinfo_);
504
505      negevals_ = wd_info_[24-1];   // Number of negative eigenvalues
506
507      if (wd_info_[0] == 0) {
508        fact_error = 0;
509      }
510      else if (wd_info_[0] == -3) {
511        /* Failure due to insufficient REAL space on a call to MA57B/BD.
512         * INFO(17) is set to a value that may suffice.  INFO(2) is set
513         * to value of LFACT.  The user can allocate a larger array and
514         * copy the contents of FACT into it using MA57E/ED, and recall
515         * MA57B/BD.
516         */
517        double  *temp;
518        ipfint  ic = 0;
519
520        wd_lfact_ = wd_info_[16];
521        temp = new double[wd_lfact_];
522
523        Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
524                       "Reallocating lfact (%d)\n", wd_lfact_);
525
526        F77_FUNC (ma57ed, MA57ED)
527        (&n, &ic, wd_keep_,
528         wd_fact_,  &wd_info_[1], temp, &wd_lfact_,
529         wd_ifact_, &wd_info_[1], NULL, &wd_lfact_,
530         wd_info_);
531
532        delete [] wd_fact_;
533        wd_fact_ = temp;
534      }
535      else if (wd_info_[0] == -4) {
536        /* Failure due to insufficient INTEGER space on a call to
537         * MA57B/BD.  INFO(18) is set to a value that may suffice.
538         * INFO(2) is set to value of LIFACT.  The user can allocate a
539         * larger array and copy the contents of IFACT into it using
540         * MA57E/ED, and recall MA57B/BD.
541         */
542
543        ipfint  *temp;
544        ipfint   ic = 1;
545
546        wd_lifact_ = wd_info_[17];
547        temp = new ipfint[wd_lifact_];
548
549        Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
550                       "Reallocating lifact (%d)\n", wd_lifact_);
551
552        F77_FUNC (ma57ed, MA57ED)
553        (&n, &ic, wd_keep_,
554         wd_fact_,  &wd_info_[1], NULL, &wd_lifact_,
555         wd_ifact_, &wd_info_[1], temp, &wd_lifact_,
556         wd_info_);
557
558        delete [] wd_ifact_;
559        wd_ifact_ = temp;
560      }
561      else if (wd_info_[0] < 0) {
562        Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
563                       "Error in MA57BD:  %d\n", wd_info_[0]);
564        Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
565                       "MA57 Error message: %s\n",
566                       ma57_err_msg[-wd_info_[1-1]].c_str());
567        return SYMSOLVER_FATAL_ERROR;
568      }
569      // Check if the system is singular.
570      else if (wd_info_[0] == 4) {
571        IpData().TimingStats().LinearSystemFactorization().End();
572        Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
573                       "System singular, rank = %d\n", wd_info_[25-1]);
574        return SYMSOLVER_SINGULAR;
575      }
576      else if (wd_info_[0] > 0) {
577        Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
578                       "Warning in MA57BD:  %d\n", wd_info_[0]);
579        Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
580                       "MA57 Warning message: %s\n",
581                       ma57_wrn_msg[wd_info_[1-1]].c_str());
582        // For now, abort the process so that we don't miss any problems
583        return SYMSOLVER_FATAL_ERROR;
584      }
585    }
586
587    // double peak_mem = 1.0e-3 * (wd_lfact*8.0 + wd_lifact*4.0 + wd_lkeep*4.0);
588
589
590    // Check whether the number of negative eigenvalues matches the
591    // requested count.
592    IpData().TimingStats().LinearSystemFactorization().End();
593    if (check_NegEVals && (numberOfNegEVals!=negevals_)) {
594      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
595                     "In Ma57TSolverInterface::Factorization: "
596                     "negevals_ = %d, but numberOfNegEVals = %d\n",
597                     negevals_, numberOfNegEVals);
598      return SYMSOLVER_WRONG_INERTIA;
599    }
600
601    return SYMSOLVER_SUCCESS;
602  }
603
604  ESymSolverStatus Ma57TSolverInterface::Backsolve(
605    Index     nrhs,
606    double    *rhs_vals)
607  {
608    DBG_START_METH("Ma27TSolverInterface::Backsolve",dbg_verbosity);
609    IpData().TimingStats().LinearSystemBackSolve().Start();
610
611    ipfint  n      = dim_;
612    ipfint  job    = 1;
613
614    ipfint  nrhs_X = nrhs;
615    ipfint  lrhs   = n;
616
617    ipfint  lwork;
618    double* work;
619
620    lwork = n * nrhs;
621    work = new double[lwork];
622
623    // For each right hand side, call MA57CD
624    // XXX MH: MA57 can do several RHSs; just do one solve...
625    // AW: Ok is the following correct?
626    if (DBG_VERBOSITY()>=2) {
627      for(Index irhs=0; irhs<nrhs; irhs++) {
628        for (Index i=0; i<dim_; i++) {
629          DBG_PRINT((2, "rhs[%2d,%5d] = %23.15e\n", irhs, i, rhs_vals[irhs*dim_+i]));
630        }
631      }
632    }
633
634    F77_FUNC (ma57cd, MA57CD)
635    (&job, &n, wd_fact_, &wd_lfact_, wd_ifact_, &wd_lifact_,
636     &nrhs_X, rhs_vals, &lrhs,
637     work, &lwork, wd_iwork_,
638     wd_icntl_, wd_info_);
639
640    if (wd_info_[0] != 0)
641      Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
642                     "Error in MA57CD:  %d.\n", wd_info_[0]);
643
644    if (DBG_VERBOSITY()>=2) {
645      for(Index irhs=0; irhs<nrhs; irhs++) {
646        for (Index i=0; i<dim_; i++) {
647          DBG_PRINT((2, "sol[%2d,%5d] = %23.15e\n", irhs, i, rhs_vals[irhs*dim_+i]));
648        }
649      }
650    }
651
652    delete [] work;
653
654    IpData().TimingStats().LinearSystemBackSolve().End();
655    return SYMSOLVER_SUCCESS;
656  }
657
658  Index Ma57TSolverInterface::NumberOfNegEVals() const
659  {
660    DBG_START_METH("Ma57TSolverInterface::NumberOfNegEVals",dbg_verbosity);
661    DBG_ASSERT(ProvidesInertia());
662    DBG_ASSERT(initialized_);
663    return negevals_;
664  }
665
666  bool Ma57TSolverInterface::IncreaseQuality()
667  {
668    DBG_START_METH("Ma57TSolverInterface::IncreaseQuality",dbg_verbosity);
669    if (pivtol_ == pivtolmax_) {
670      return false;
671    }
672    pivtol_changed_ = true;
673
674    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
675                   "Indreasing pivot tolerance for MA57 from %7.2e ",
676                   pivtol_);
677    pivtol_ = Min(pivtolmax_, pow(pivtol_,0.75));
678    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
679                   "to %7.2e.\n",
680                   pivtol_);
681    return true;
682  }
683
684} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.