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

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

added dummy ipoma.f file in case CUTEr is not installed

  • 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 568 2005-11-01 23:28:48Z 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      wd_keep_(NULL),
241      wd_iwork_(NULL),
242      wd_fact_(NULL),
243      wd_ifact_(NULL),
244      a_(NULL)
245  {
246    DBG_START_METH("Ma57TSolverInterface::Ma57TSolverInterface()",
247                   dbg_verbosity);
248  }
249
250  Ma57TSolverInterface::~Ma57TSolverInterface()
251  {
252    DBG_START_METH("Ma57TSolverInterface::~Ma57TSolverInterface()",
253                   dbg_verbosity);
254    delete [] a_;
255
256    delete [] wd_fact_;
257    delete [] wd_ifact_;
258
259    delete [] wd_iwork_;
260    delete [] wd_keep_;
261  }
262
263  void Ma57TSolverInterface::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
264  {
265    roptions->AddBoundedNumberOption(
266      "ma57_pivtol",
267      "Pivot tolerance for the linear solver MA57.",
268      0.0, true, 1.0, true, 1e-8,
269      "A smaller number pivots for sparsity, "
270      "a larger number pivots for stability.");
271    roptions->AddBoundedNumberOption(
272      "ma57_pivtolmax",
273      "Maximum pivot tolerance for the linear solver MA57.",
274      0.0, true, 1.0, true, 1e-4,
275      "Ipopt may increase pivtol as high as pivtolmax "
276      "to get a more accurate solution to the linear system.");
277    roptions->AddLowerBoundedNumberOption(
278      "ma57_pre_alloc",
279      "Safety factor for work space memory allocation for the linear solver MA57.",
280      1., false, 3.,
281      "If 1 is chosen, the suggested amount of work space is used.  However, "
282      "choosen a larger number might avoid reallocation if the suggest values "
283      "do not suffice.");
284  }
285
286  bool Ma57TSolverInterface::InitializeImpl(const OptionsList&  options,
287      const std::string&    prefix)
288  {
289    // Obtain the options settings
290    options.GetNumericValue("ma57_pivtol", pivtol_, prefix);
291    if(options.GetNumericValue("ma57_pivtolmax", pivtolmax_, prefix)) {
292      ASSERT_EXCEPTION(pivtolmax_>=pivtol_, OPTION_INVALID,
293                       "Option \"pivtolmax\": This value must be between "
294                       "pivtol and 1.");
295    }
296    else {
297      pivtolmax_ = Max(pivtolmax_, pivtol_);
298    }
299
300    options.GetNumericValue("ma57_pre_alloc", ma57_pre_alloc_, prefix);
301
302    // The following option is registered by OrigIpoptNLP
303    options.GetBoolValue("warm_start_same_structure",
304                         warm_start_same_structure_, prefix);
305    DBG_ASSERT(!warm_start_same_structure_ && "warm_start_same_structure not yet implemented");
306
307
308    /* Initialize. */
309    F77_FUNC (ma57id, MA57ID) (wd_cntl_, wd_icntl_);
310
311    /* Custom settings for MA57. */
312    wd_icntl_[1-1] = 0;      /* Error stream */
313    wd_icntl_[2-1] = 0;      /* Warning stream. */
314
315    wd_icntl_[4-1] = 1;      /* Print statistics.  NOT Used. */
316    wd_icntl_[5-1] = 0;      /* Print error. */
317
318    // wd_icntl[6-1] = 0;       /* Pivoting order. */
319
320    wd_cntl_[1-1]  = pivtol_;    /* Pivot threshold. */
321    wd_icntl_[7-1] = 1;      /* Pivoting strategy. */
322
323    // wd_icntl[8-1] = 0;       /* Retry factorization. */
324
325    if (!warm_start_same_structure_) {
326      dim_=0;
327      nonzeros_=0;
328      delete [] a_;
329      delete [] wd_fact_;
330      delete [] wd_ifact_;
331      delete [] wd_iwork_;
332      delete [] wd_keep_;
333    }
334    else {
335      ASSERT_EXCEPTION(dim_>0 && nonzeros_>0, INVALID_WARMSTART,
336                       "Ma57TSolverInterface called with warm_start_same_structure, "
337                       "but the problem is solved for the first time.");
338    }
339
340    return true;
341  }
342
343  ESymSolverStatus
344  Ma57TSolverInterface::MultiSolve(bool         new_matrix,
345                                   const Index*     airn,
346                                   const Index*     ajcn,
347                                   Index        nrhs,
348                                   double*      rhs_vals,
349                                   bool         check_NegEVals,
350                                   Index        numberOfNegEVals)
351  {
352    DBG_START_METH("Ma57TSolverInterface::MultiSolve",dbg_verbosity);
353
354    // DBG_ASSERT(!check_NegEVals || ProvidesInertia());
355    // DBG_ASSERT(initialized_);
356    // DBG_ASSERT(la_!=0);
357
358    if (pivtol_changed_) {
359      DBG_PRINT((1,"Pivot tolerance has changed.\n"));
360      pivtol_changed_ = false;
361      // If the pivot tolerance has been changed but the matrix is not
362      // new, we have to request the values for the matrix again to do
363      // the factorization again.
364      if (!new_matrix) {
365        DBG_PRINT((1,"Ask caller to call again.\n"));
366        refactorize_ = true;
367        return SYMSOLVER_CALL_AGAIN;
368      }
369    }
370
371    // check if a factorization has to be done
372    DBG_PRINT((1, "new_matrix = %d\n", new_matrix));
373    if (new_matrix || refactorize_) {
374      // perform the factorization
375      ESymSolverStatus retval;
376      retval = Factorization(airn, ajcn, check_NegEVals, numberOfNegEVals);
377      if (retval!=SYMSOLVER_SUCCESS) {
378        DBG_PRINT((1, "FACTORIZATION FAILED!\n"));
379        return retval;  // Matrix singular or error occurred
380      }
381      refactorize_ = false;
382    }
383
384    // do the backsolve
385    return Backsolve(nrhs, rhs_vals);
386  }
387
388  double* Ma57TSolverInterface::GetValuesArrayPtr()
389  {
390    DBG_START_METH("Ma57TSolverInterface::GetValuesArrayPtr",dbg_verbosity);
391    DBG_ASSERT(initialized_);
392
393    return a_;
394  }
395
396  /** Initialize the local copy of the positions of the nonzero
397      elements */
398  ESymSolverStatus
399  Ma57TSolverInterface::InitializeStructure(
400    Index       dim,
401    Index     nonzeros,
402    const Index*  airn,
403    const Index*  ajcn)
404  {
405    DBG_START_METH("Ma57TSolverInterface::InitializeStructure",dbg_verbosity);
406
407    ESymSolverStatus retval = SYMSOLVER_SUCCESS;
408    if (!warm_start_same_structure_) {
409      dim_ = dim;
410      nonzeros_ = nonzeros;
411      // for MA57, a_ only has to be as long as the number of nonzero
412      // elements
413      delete [] a_;
414      a_ = new double [nonzeros_];
415
416      // Do the symbolic facotrization
417      retval = SymbolicFactorization(airn, ajcn);
418      if (retval != SYMSOLVER_SUCCESS ) {
419        return retval;
420      }
421    }
422    else {
423      ASSERT_EXCEPTION(dim_==dim && nonzeros_==nonzeros, INVALID_WARMSTART,
424                       "Ma57TSolverInterface called with warm_start_same_structure, "
425                       "but the problem size has changed.");
426    }
427
428    initialized_ = true;
429
430    return retval;
431  }
432
433  ESymSolverStatus
434  Ma57TSolverInterface::SymbolicFactorization(const Index*  airn,
435      const Index*  ajcn)
436  {
437    DBG_START_METH("Ma57TSolverInterface::SymbolicFactorization",dbg_verbosity);
438    IpData().TimingStats().LinearSystemSymbolicFactorization().Start();
439
440    ipfint n  = dim_;
441    ipfint ne = nonzeros_;
442
443    wd_lkeep_ = 5*n + ne + Max (n,ne) + 42;
444
445    wd_iwork_ = new ipfint[5*n];
446    wd_keep_  = new ipfint[wd_lkeep_];
447
448    F77_FUNC (ma57ad, MA57AD)
449    (&n, &ne, airn, ajcn, &wd_lkeep_, wd_keep_, wd_iwork_,
450     wd_icntl_, wd_info_, wd_rinfo_);
451
452    if (wd_info_[0] < 0) {
453      Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
454                     "*** Error from MA57AD *** INFO(0) = %d\n", wd_info_[0]);
455    }
456
457    wd_lfact_  = (ipfint)((Number)wd_info_[8] * ma57_pre_alloc_);
458    wd_lifact_ = (ipfint)((Number)wd_info_[9] * ma57_pre_alloc_);
459
460    // XXX MH:  Why is this necessary?  Is `::Factorization' called more
461    // than once per object lifetime?  Where should allocation take
462    // place, then?
463
464    // AW: I moved this here now - my understanding is that wd_info[8]
465    // and wd_info[9] are computed here, so we can just allocate the
466    // amount of memory here.  I don't think there is any need to
467    // reallocate it later for every factorization
468    delete [] wd_fact_;
469    delete [] wd_ifact_;
470
471    wd_fact_  = new double[wd_lfact_];
472    wd_ifact_ = new int[wd_lifact_];
473
474    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
475                   "Suggested lfact  (*%e):  %d\n", ma57_pre_alloc_, wd_lfact_);
476    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
477                   "Suggested lifact (*%e):  %d\n", ma57_pre_alloc_, wd_lifact_);
478
479    IpData().TimingStats().LinearSystemSymbolicFactorization().End();
480    return SYMSOLVER_SUCCESS;
481  }
482
483  ESymSolverStatus
484  Ma57TSolverInterface::Factorization(const Index*  airn,
485                                      const Index*  ajcn,
486                                      bool          check_NegEVals,
487                                      Index         numberOfNegEVals)
488  {
489    DBG_START_METH("Ma57TSolverInterface::Factorization",dbg_verbosity);
490    IpData().TimingStats().LinearSystemFactorization().Start();
491
492    int fact_error = 1;
493
494    ipfint n  = dim_;
495    ipfint ne = nonzeros_;
496
497    while (fact_error > 0) {
498      F77_FUNC (ma57bd, MA57BD)
499      (&n, &ne, a_, wd_fact_, &wd_lfact_, wd_ifact_, &wd_lifact_,
500       &wd_lkeep_, wd_keep_, wd_iwork_,
501       wd_icntl_, wd_cntl_, wd_info_, wd_rinfo_);
502
503      negevals_ = wd_info_[24-1];   // Number of negative eigenvalues
504
505      if (wd_info_[0] == 0) {
506        fact_error = 0;
507      }
508      else if (wd_info_[0] == -3) {
509        /* Failure due to insufficient REAL space on a call to MA57B/BD.
510         * INFO(17) is set to a value that may suffice.  INFO(2) is set
511         * to value of LFACT.  The user can allocate a larger array and
512         * copy the contents of FACT into it using MA57E/ED, and recall
513         * MA57B/BD.
514         */
515        double  *temp;
516        ipfint  ic = 0;
517
518        wd_lfact_ = wd_info_[16];
519        temp = new double[wd_lfact_];
520
521        Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
522                       "Reallocating lfact (%d)\n", wd_lfact_);
523
524        F77_FUNC (ma57ed, MA57ED)
525        (&n, &ic, wd_keep_,
526         wd_fact_,  &wd_info_[1], temp, &wd_lfact_,
527         wd_ifact_, &wd_info_[1], NULL, &wd_lfact_,
528         wd_info_);
529
530        delete [] wd_fact_;
531        wd_fact_ = temp;
532      }
533      else if (wd_info_[0] == -4) {
534        /* Failure due to insufficient INTEGER space on a call to
535         * MA57B/BD.  INFO(18) is set to a value that may suffice.
536         * INFO(2) is set to value of LIFACT.  The user can allocate a
537         * larger array and copy the contents of IFACT into it using
538         * MA57E/ED, and recall MA57B/BD.
539         */
540
541        ipfint  *temp;
542        ipfint   ic = 1;
543
544        wd_lifact_ = wd_info_[17];
545        temp = new ipfint[wd_lifact_];
546
547        Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
548                       "Reallocating lifact (%d)\n", wd_lifact_);
549
550        F77_FUNC (ma57ed, MA57ED)
551        (&n, &ic, wd_keep_,
552         wd_fact_,  &wd_info_[1], NULL, &wd_lifact_,
553         wd_ifact_, &wd_info_[1], temp, &wd_lifact_,
554         wd_info_);
555
556        delete [] wd_ifact_;
557        wd_ifact_ = temp;
558      }
559      else if (wd_info_[0] < 0) {
560        Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
561                       "Error in MA57BD:  %d\n", wd_info_[0]);
562        Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
563                       "MA57 Error message: %s\n",
564                       ma57_err_msg[-wd_info_[1-1]].c_str());
565        return SYMSOLVER_FATAL_ERROR;
566      }
567      // Check if the system is singular.
568      else if (wd_info_[0] == 4) {
569        IpData().TimingStats().LinearSystemFactorization().End();
570        Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
571                       "System singular, rank = %d\n", wd_info_[25-1]);
572        return SYMSOLVER_SINGULAR;
573      }
574      else if (wd_info_[0] > 0) {
575        Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
576                       "Warning in MA57BD:  %d\n", wd_info_[0]);
577        Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
578                       "MA57 Warning message: %s\n",
579                       ma57_wrn_msg[wd_info_[1-1]].c_str());
580        // For now, abort the process so that we don't miss any problems
581        return SYMSOLVER_FATAL_ERROR;
582      }
583    }
584
585    // double peak_mem = 1.0e-3 * (wd_lfact*8.0 + wd_lifact*4.0 + wd_lkeep*4.0);
586
587
588    // Check whether the number of negative eigenvalues matches the
589    // requested count.
590    IpData().TimingStats().LinearSystemFactorization().End();
591    if (check_NegEVals && (numberOfNegEVals!=negevals_)) {
592      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
593                     "In Ma57TSolverInterface::Factorization: "
594                     "negevals_ = %d, but numberOfNegEVals = %d\n",
595                     negevals_, numberOfNegEVals);
596      return SYMSOLVER_WRONG_INERTIA;
597    }
598
599    return SYMSOLVER_SUCCESS;
600  }
601
602  ESymSolverStatus Ma57TSolverInterface::Backsolve(
603    Index     nrhs,
604    double    *rhs_vals)
605  {
606    DBG_START_METH("Ma27TSolverInterface::Backsolve",dbg_verbosity);
607    IpData().TimingStats().LinearSystemBackSolve().Start();
608
609    ipfint  n      = dim_;
610    ipfint  job    = 1;
611
612    ipfint  nrhs_X = nrhs;
613    ipfint  lrhs   = n;
614
615    ipfint  lwork;
616    double* work;
617
618    lwork = n * nrhs;
619    work = new double[lwork];
620
621    // For each right hand side, call MA57CD
622    // XXX MH: MA57 can do several RHSs; just do one solve...
623    // AW: Ok is the following correct?
624    if (DBG_VERBOSITY()>=2) {
625      for(Index irhs=0; irhs<nrhs; irhs++) {
626        for (Index i=0; i<dim_; i++) {
627          DBG_PRINT((2, "rhs[%2d,%5d] = %23.15e\n", irhs, i, rhs_vals[irhs*dim_+i]));
628        }
629      }
630    }
631
632    F77_FUNC (ma57cd, MA57CD)
633    (&job, &n, wd_fact_, &wd_lfact_, wd_ifact_, &wd_lifact_,
634     &nrhs_X, rhs_vals, &lrhs,
635     work, &lwork, wd_iwork_,
636     wd_icntl_, wd_info_);
637
638    if (wd_info_[0] != 0)
639      Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
640                     "Error in MA57CD:  %d.\n", wd_info_[0]);
641
642    if (DBG_VERBOSITY()>=2) {
643      for(Index irhs=0; irhs<nrhs; irhs++) {
644        for (Index i=0; i<dim_; i++) {
645          DBG_PRINT((2, "sol[%2d,%5d] = %23.15e\n", irhs, i, rhs_vals[irhs*dim_+i]));
646        }
647      }
648    }
649
650    delete [] work;
651
652    IpData().TimingStats().LinearSystemBackSolve().End();
653    return SYMSOLVER_SUCCESS;
654  }
655
656  Index Ma57TSolverInterface::NumberOfNegEVals() const
657  {
658    DBG_START_METH("Ma57TSolverInterface::NumberOfNegEVals",dbg_verbosity);
659    DBG_ASSERT(ProvidesInertia());
660    DBG_ASSERT(initialized_);
661    return negevals_;
662  }
663
664  bool Ma57TSolverInterface::IncreaseQuality()
665  {
666    DBG_START_METH("Ma57TSolverInterface::IncreaseQuality",dbg_verbosity);
667    if (pivtol_ == pivtolmax_) {
668      return false;
669    }
670    pivtol_changed_ = true;
671
672    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
673                   "Indreasing pivot tolerance for MA57 from %7.2e ",
674                   pivtol_);
675    pivtol_ = Min(pivtolmax_, pow(pivtol_,0.75));
676    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
677                   "to %7.2e.\n",
678                   pivtol_);
679    return true;
680  }
681
682} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.