source: branches/dev/Algorithm/IpPDFullSpaceSolver.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: 29.5 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: IpPDFullSpaceSolver.cpp 529 2005-09-29 21:12:38Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8
9#include "IpPDFullSpaceSolver.hpp"
10#include "IpDebug.hpp"
11
12namespace Ipopt
13{
14
15#ifdef IP_DEBUG
16  static const Index dbg_verbosity = 0;
17#endif
18
19  PDFullSpaceSolver::PDFullSpaceSolver(AugSystemSolver& augSysSolver,
20                                       PDPerturbationHandler& perturbHandler)
21      :
22      PDSystemSolver(),
23      augSysSolver_(&augSysSolver),
24      perturbHandler_(&perturbHandler),
25      dummy_cache_(1)
26  {
27    DBG_START_METH("PDFullSpaceSolver::PDFullSpaceSolver",dbg_verbosity);
28  }
29
30  PDFullSpaceSolver::~PDFullSpaceSolver()
31  {
32    DBG_START_METH("PDFullSpaceSolver::~PDFullSpaceSolver()",dbg_verbosity);
33  }
34
35  void PDFullSpaceSolver::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
36  {
37    roptions->AddLowerBoundedIntegerOption(
38      "min_refinement_steps",
39      "Minimum number of iterative refinement steps per linear system solve.",
40      0, 1,
41      "Iterative refinement (on the full unsymmetric system) is performed for "
42      "each right hand side.  This option determines the minimum number "
43      "of iterative refinements (i.e. at least \"min_refinement_steps\" "
44      "iterative refinement steps are enforced per right hand side.)");
45    roptions->AddLowerBoundedIntegerOption(
46      "max_refinement_steps",
47      "Maximum number of iterative refinement steps per linear system solve.",
48      0, 10,
49      "Iterative refinement (on the full unsymmetric system) is performed for "
50      "each right hand side.  This option determines the maximum number "
51      "of iterative refinement steps.");
52    roptions->AddLowerBoundedNumberOption(
53      "residual_ratio_max",
54      "Iterative refinement tolerance",
55      0.0, true, 1e-10,
56      "Iterative refinement is performed until the residual test ratio is "
57      "less than this tolerance (or until \"max_refinement_steps\" refinement "
58      "steps are performed).");
59    roptions->AddLowerBoundedNumberOption(
60      "residual_ratio_singular",
61      "Threshold for declaring linear system singular after failed iterative refinement.",
62      0.0, true, 1e-5,
63      "If the residual test ratio is larger than this value after failed "
64      "iterative refinement, the algorithm pretends that the linear system is "
65      "singular.");
66    // ToDo Think about following option - are the correct norms used?
67    roptions->AddLowerBoundedNumberOption(
68      "residual_improvement_factor",
69      "Minimal required reduction of residual test ratio in iterative refinement.",
70      0.0, true, 0.999999999,
71      "If the improvement of the residual test ratio made by one iterative "
72      "refinement step is not better than this factor, iterative refinement "
73      "is aborted.");
74  }
75
76
77  bool PDFullSpaceSolver::InitializeImpl(const OptionsList& options,
78                                         const std::string& prefix)
79  {
80    // Check for the algorithm options
81    options.GetIntegerValue("min_refinement_steps", min_refinement_steps_, prefix);
82    options.GetIntegerValue("max_refinement_steps", max_refinement_steps_, prefix);
83    ASSERT_EXCEPTION(max_refinement_steps_ >= min_refinement_steps_, OPTION_INVALID,
84                     "Option \"max_refinement_steps\": This value must be larger than or equal to min_refinement_steps (default 1)");
85
86    options.GetNumericValue("residual_ratio_max", residual_ratio_max_, prefix);
87    options.GetNumericValue("residual_ratio_singular", residual_ratio_singular_, prefix);
88    ASSERT_EXCEPTION(residual_ratio_singular_ > residual_ratio_max_, OPTION_INVALID,
89                     "Option \"residual_ratio_singular\": This value must be larger than residual_ratio_max.");
90    options.GetNumericValue("residual_improvement_factor", residual_improvement_factor_, prefix);
91
92    // Reset internal flags and data
93    augsys_improved_ = false;
94
95    if (!augSysSolver_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
96                                   options, prefix)) {
97      return false;
98    }
99
100    return perturbHandler_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
101                                       options, prefix);
102  }
103
104  void PDFullSpaceSolver::Solve(Number alpha,
105                                Number beta,
106                                const IteratesVector& rhs,
107                                IteratesVector& res,
108                                bool allow_inexact,
109                                bool improve_solution /* = false */)
110  {
111    DBG_START_METH("PDFullSpaceSolver::Solve",dbg_verbosity);
112    DBG_ASSERT(!allow_inexact || !improve_solution);
113    DBG_ASSERT(!improve_solution || beta==0.);
114
115    // Timing of PDSystem solver starts here
116    IpData().TimingStats().PDSystemSolverTotal().Start();
117
118    DBG_PRINT_VECTOR(2, "rhs_x", *rhs.x());
119    DBG_PRINT_VECTOR(2, "rhs_s", *rhs.s());
120    DBG_PRINT_VECTOR(2, "rhs_c", *rhs.y_c());
121    DBG_PRINT_VECTOR(2, "rhs_d", *rhs.y_d());
122    DBG_PRINT_VECTOR(2, "rhs_zL", *rhs.z_L());
123    DBG_PRINT_VECTOR(2, "rhs_zU", *rhs.z_U());
124    DBG_PRINT_VECTOR(2, "rhs_vL", *rhs.v_L());
125    DBG_PRINT_VECTOR(2, "rhs_vU", *rhs.v_U());
126    DBG_PRINT_VECTOR(2, "res_x in", *res.x());
127    DBG_PRINT_VECTOR(2, "res_s in", *res.s());
128    DBG_PRINT_VECTOR(2, "res_c in", *res.y_c());
129    DBG_PRINT_VECTOR(2, "res_d in", *res.y_d());
130    DBG_PRINT_VECTOR(2, "res_zL in", *res.z_L());
131    DBG_PRINT_VECTOR(2, "res_zU in", *res.z_U());
132    DBG_PRINT_VECTOR(2, "res_vL in", *res.v_L());
133    DBG_PRINT_VECTOR(2, "res_vU in", *res.v_U());
134
135    // if beta is nonzero, keep a copy of the incoming values in res_ */
136    SmartPtr<IteratesVector> copy_res;
137    if (beta != 0.) {
138      copy_res = res.MakeNewIteratesVectorCopy();
139    }
140
141    // Receive data about matrix
142    SmartPtr<const Vector> x = IpData().curr()->x();
143    SmartPtr<const Vector> s = IpData().curr()->s();
144    SmartPtr<const SymMatrix> W = IpData().W();
145    SmartPtr<const Matrix> J_c = IpCq().curr_jac_c();
146    SmartPtr<const Matrix> J_d = IpCq().curr_jac_d();
147    SmartPtr<const Matrix> Px_L = IpNLP().Px_L();
148    SmartPtr<const Matrix> Px_U = IpNLP().Px_U();
149    SmartPtr<const Matrix> Pd_L = IpNLP().Pd_L();
150    SmartPtr<const Matrix> Pd_U = IpNLP().Pd_U();
151    SmartPtr<const Vector> z_L = IpData().curr()->z_L();
152    SmartPtr<const Vector> z_U = IpData().curr()->z_U();
153    SmartPtr<const Vector> v_L = IpData().curr()->v_L();
154    SmartPtr<const Vector> v_U = IpData().curr()->v_U();
155    SmartPtr<const Vector> slack_x_L = IpCq().curr_slack_x_L();
156    SmartPtr<const Vector> slack_x_U = IpCq().curr_slack_x_U();
157    SmartPtr<const Vector> slack_s_L = IpCq().curr_slack_s_L();
158    SmartPtr<const Vector> slack_s_U = IpCq().curr_slack_s_U();
159    SmartPtr<const Vector> sigma_x = IpCq().curr_sigma_x();
160    SmartPtr<const Vector> sigma_s = IpCq().curr_sigma_s();
161    DBG_PRINT_VECTOR(2, "Sigma_x", *sigma_x);
162    DBG_PRINT_VECTOR(2, "Sigma_s", *sigma_s);
163
164    bool done = false;
165    // The following flag is set to true, if we asked the linear
166    // solver to improve the quality of the solution in
167    // the next solve
168    bool resolve_with_better_quality = false;
169    // the following flag is set to true, if iterative refinement
170    // failed and we want to try if a modified system is able to
171    // remedy that problem by pretending the matrix is singular
172    bool pretend_singular = false;
173    bool pretend_singular_last_time = false;
174
175    // Beginning of loop for solving the system (including all
176    // modifications for the linear system to ensure good solution
177    // quality)
178    while (!done) {
179
180      // if improve_solution is true, we are given already a solution
181      // from the calling function, so we can skip the first solve
182      bool solve_retval = true;
183      if (!improve_solution) {
184        solve_retval =
185          SolveOnce(resolve_with_better_quality, pretend_singular,
186                    *W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U,
187                    *v_L, *v_U, *slack_x_L, *slack_x_U, *slack_s_L, *slack_s_U,
188                    *sigma_x, *sigma_s, 1., 0., rhs, res);
189        resolve_with_better_quality = false;
190        pretend_singular = false;
191      }
192      improve_solution = false;
193
194      if (!solve_retval) {
195        // If system seems not to be solvable, we set the search
196        // direction to zero, and hope that the line search will take
197        // care of this (e.g. call the restoration phase).  ToDo: We
198        // might want to use a more explicit cue later.
199        res.Set(0.0);
200        IpData().TimingStats().PDSystemSolverTotal().End();
201        return;
202      }
203
204      if (allow_inexact) {
205        // no safety checks required
206        break;
207      }
208
209      // Get space for the residual
210      SmartPtr<IteratesVector> resid = res.MakeNewIteratesVector(true);
211
212      // ToDo don't to that after max refinement?
213      ComputeResiduals(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U,
214                       *z_L, *z_U, *v_L, *v_U, *slack_x_L, *slack_x_U,
215                       *slack_s_L, *slack_s_U, *sigma_x, *sigma_s,
216                       alpha, beta, rhs, res, *resid);
217
218      Number residual_ratio =
219        ComputeResidualRatio(rhs, res, *resid);
220      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
221                     "residual_ratio = %e\n", residual_ratio);
222      Number residual_ratio_old = residual_ratio;
223
224      // Beginning of loop for iterative refinement
225      Index num_iter_ref = 0;
226      bool quit_refinement = false;
227      while (!allow_inexact && !quit_refinement &&
228             (num_iter_ref < min_refinement_steps_ ||
229              residual_ratio > residual_ratio_max_) ) {
230
231        // To the next back solve
232        solve_retval =
233          SolveOnce(resolve_with_better_quality, false,
234                    *W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U,
235                    *v_L, *v_U, *slack_x_L, *slack_x_U, *slack_s_L, *slack_s_U,
236                    *sigma_x, *sigma_s, -1., 1., *resid, res);
237        ASSERT_EXCEPTION(solve_retval, INTERNAL_ABORT,
238                         "SolveOnce returns false.");
239
240        ComputeResiduals(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U,
241                         *z_L, *z_U, *v_L, *v_U, *slack_x_L, *slack_x_U,
242                         *slack_s_L, *slack_s_U, *sigma_x, *sigma_s,
243                         alpha, beta, rhs, res, *resid);
244
245        residual_ratio =
246          ComputeResidualRatio(rhs, res, *resid);
247        Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
248                       "residual_ratio = %e\n", residual_ratio);
249
250        num_iter_ref++;
251        // Check if we have to give up on iterative refinement
252        if (num_iter_ref>min_refinement_steps_ &&
253            (num_iter_ref>max_refinement_steps_ ||
254             residual_ratio>residual_improvement_factor_*residual_ratio_old)) {
255
256          Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
257                         "Iterative refinement failed with residual_ratio = %e\n", residual_ratio);
258          quit_refinement = true;
259
260          // Pretend singularity only once - if it didn't help, we
261          // have to live with what we got so far
262          resolve_with_better_quality = false;
263          DBG_PRINT((1, "pretend_singular = %d\n", pretend_singular));
264          if (!pretend_singular_last_time) {
265            // First try if we can ask the augmented system solver to
266            // improve the quality of the solution (only if that hasn't
267            // been done before for this linear system)
268            if (!augsys_improved_) {
269              Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
270                             "Asking augmented system solver to improve quality of its solutions.\n");
271              augsys_improved_ = augSysSolver_->IncreaseQuality();
272              if (augsys_improved_) {
273                IpData().Append_info_string("q");
274                resolve_with_better_quality = true;
275              }
276              else {
277                // solver said it cannot improve quality, so let
278                // possibly conclude that the current modification is
279                // singular
280                pretend_singular = true;
281              }
282            }
283            else {
284              // we had already asked the solver before to improve the
285              // quality of the solution, so let's now pretend that the
286              // modification is possibly singular
287              pretend_singular = true;
288            }
289            pretend_singular_last_time = pretend_singular;
290            if (pretend_singular) {
291              // let's only conclude that the current linear system
292              // including modifications is singular, if the residual is
293              // quite bad
294              if (residual_ratio < residual_ratio_singular_) {
295                pretend_singular = false;
296                IpData().Append_info_string("S");
297                Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
298                               "Just accept current solution.\n");
299              }
300              else {
301                IpData().Append_info_string("s");
302                Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
303                               "Pretend that the current system (including modifications) is singular.\n");
304              }
305            }
306          }
307          else {
308            pretend_singular = false;
309            DBG_PRINT((1,"Resetting pretend_singular to false.\n"));
310          }
311        }
312
313        residual_ratio_old = residual_ratio;
314      } // End of loop for iterative refinement
315
316      done = !(resolve_with_better_quality) && !(pretend_singular);
317
318    } // End of loop for solving the linear system (incl. modifications)
319
320    // Finally let's assemble the res result vectors
321    if (alpha != 0.) {
322      res.Scal(alpha);
323    }
324
325    if (beta != 0.) {
326      res.Axpy(beta, *copy_res);
327    }
328
329    DBG_PRINT_VECTOR(2, "res_x", *res.x());
330    DBG_PRINT_VECTOR(2, "res_s", *res.s());
331    DBG_PRINT_VECTOR(2, "res_c", *res.y_c());
332    DBG_PRINT_VECTOR(2, "res_d", *res.y_d());
333    DBG_PRINT_VECTOR(2, "res_zL", *res.z_L());
334    DBG_PRINT_VECTOR(2, "res_zU", *res.z_U());
335    DBG_PRINT_VECTOR(2, "res_vL", *res.v_L());
336    DBG_PRINT_VECTOR(2, "res_vU", *res.v_U());
337
338    IpData().TimingStats().PDSystemSolverTotal().End();
339  }
340
341  bool PDFullSpaceSolver::SolveOnce(bool resolve_with_better_quality,
342                                    bool pretend_singular,
343                                    const SymMatrix& W,
344                                    const Matrix& J_c,
345                                    const Matrix& J_d,
346                                    const Matrix& Px_L,
347                                    const Matrix& Px_U,
348                                    const Matrix& Pd_L,
349                                    const Matrix& Pd_U,
350                                    const Vector& z_L,
351                                    const Vector& z_U,
352                                    const Vector& v_L,
353                                    const Vector& v_U,
354                                    const Vector& slack_x_L,
355                                    const Vector& slack_x_U,
356                                    const Vector& slack_s_L,
357                                    const Vector& slack_s_U,
358                                    const Vector& sigma_x,
359                                    const Vector& sigma_s,
360                                    Number alpha,
361                                    Number beta,
362                                    const IteratesVector& rhs,
363                                    IteratesVector& res)
364  {
365    // TO DO LIST:
366    //
367    // 1. decide for reasonable return codes (e.g. fatal error, too
368    //    ill-conditioned...)
369    // 2. Make constants parameters that can be set from the outside
370    // 3. Get Information out of Ipopt structures
371    // 4. add heuristic for structurally singular problems
372    // 5. see if it makes sense to distinguish delta_x and delta_s,
373    //    or delta_c and delta_d
374    // 6. increase pivot tolerance if number of get evals so too small
375    DBG_START_METH("PDFullSpaceSolver::SolveOnce",dbg_verbosity);
376
377    IpData().TimingStats().PDSystemSolverSolveOnce().Start();
378
379    // Compute the right hand side for the augmented system formulation
380    SmartPtr<Vector> augRhs_x = rhs.x()->MakeNewCopy();
381    Px_L.AddMSinvZ(1.0, slack_x_L, *rhs.z_L(), *augRhs_x);
382    Px_U.AddMSinvZ(-1.0, slack_x_U, *rhs.z_U(), *augRhs_x);
383
384    SmartPtr<Vector> augRhs_s = rhs.s()->MakeNewCopy();
385    Pd_L.AddMSinvZ(1.0, slack_s_L, *rhs.v_L(), *augRhs_s);
386    Pd_U.AddMSinvZ(-1.0, slack_s_U, *rhs.v_U(), *augRhs_s);
387
388    // Get space into which we can put the solution of the augmented system
389    SmartPtr<IteratesVector> sol = res.MakeNewIteratesVector(true);
390
391    // Now check whether any data has changed
392    std::vector<const TaggedObject*> deps(13);
393    deps[0] = &W;
394    deps[1] = &J_c;
395    deps[2] = &J_d;
396    deps[3] = &z_L;
397    deps[4] = &z_U;
398    deps[5] = &v_L;
399    deps[6] = &v_U;
400    deps[7] = &slack_x_L;
401    deps[8] = &slack_x_U;
402    deps[9] = &slack_s_L;
403    deps[10] = &slack_s_U;
404    deps[11] = &sigma_x;
405    deps[12] = &sigma_s;
406    void* dummy;
407    bool uptodate = dummy_cache_.GetCachedResult(dummy, deps);
408    if (!uptodate) {
409      dummy_cache_.AddCachedResult(dummy, deps);
410      augsys_improved_ = false;
411    }
412    // improve_current_solution can only be true, if that system has
413    // been solved before
414    DBG_ASSERT((!resolve_with_better_quality && !pretend_singular) || uptodate);
415
416    ESymSolverStatus retval;
417    if (uptodate && !pretend_singular) {
418
419      // Get the perturbation values
420      Number delta_x;
421      Number delta_s;
422      Number delta_c;
423      Number delta_d;
424      perturbHandler_->CurrentPerturbation(delta_x, delta_s, delta_c, delta_d);
425
426      // No need to go throught the pain of finding the appropriate
427      // values for the deltas, because the matrix hasn't changed since
428      // the last call.  So, just call the Solve Method
429      //
430      // Note: resolve_with_better_quality is true, then the Solve
431      // method has already asked the augSysSolver to increase the
432      // quality at the end solve, and we are now getting the solution
433      // with that better quality
434      retval = augSysSolver_->Solve(&W, &sigma_x, delta_x,
435                                    &sigma_s, delta_s, &J_c, NULL,
436                                    delta_c, &J_d, NULL, delta_d,
437                                    *augRhs_x, *augRhs_s, *rhs.y_c(), *rhs.y_d(),
438                                    *sol->x_NonConst(), *sol->s_NonConst(),
439                                    *sol->y_c_NonConst(), *sol->y_d_NonConst(),
440                                    false, 0);
441      if (retval!=SYMSOLVER_SUCCESS) {
442        IpData().TimingStats().PDSystemSolverSolveOnce().End();
443        return false;
444      }
445    }
446    else {
447      Index numberOfEVals=rhs.y_c()->Dim()+rhs.y_d()->Dim();
448      // counter for the number of trial evaluations
449      // (ToDo is not at the correct place)
450      Index count = 0;
451
452      // Get the very first perturbation values from the perturbation
453      // Handler
454      Number delta_x;
455      Number delta_s;
456      Number delta_c;
457      Number delta_d;
458      perturbHandler_->ConsiderNewSystem(delta_x, delta_s, delta_c, delta_d);
459
460      retval = SYMSOLVER_SINGULAR;
461      bool fail = false;
462
463      while (retval!= SYMSOLVER_SUCCESS && !fail) {
464
465        if (pretend_singular) {
466          retval = SYMSOLVER_SINGULAR;
467          pretend_singular = false;
468        }
469        else {
470          count++;
471          Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
472                         "Solving system with delta_x=%e delta_s=%e\n                    delta_c=%e delta_d=%e\n",
473                         delta_x, delta_s, delta_c, delta_d);
474          retval = augSysSolver_->Solve(&W, &sigma_x, delta_x,
475                                        &sigma_s, delta_s, &J_c, NULL,
476                                        delta_c, &J_d, NULL, delta_d,
477                                        *augRhs_x, *augRhs_s, *rhs.y_c(), *rhs.y_d(),
478                                        *sol->x_NonConst(), *sol->s_NonConst(),
479                                        *sol->y_c_NonConst(), *sol->y_d_NonConst(),
480                                        true, numberOfEVals);
481        }
482        assert(retval!=SYMSOLVER_FATAL_ERROR); //TODO make return code
483        if (retval==SYMSOLVER_SINGULAR &&
484            (rhs.y_c()->Dim()+rhs.y_d()->Dim() > 0) ) {
485
486          // Get new perturbation factors from the perturbation
487          // handlers for the singular case
488          bool pert_return =
489                             perturbHandler_->PerturbForSingularity(delta_x, delta_s,
490                                                                    delta_c, delta_d);
491          if (!pert_return) {
492            Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
493                           "PerturbForSingularity can't be done\n");
494            IpData().TimingStats().PDSystemSolverSolveOnce().End();
495            return false;
496          }
497        }
498        else if (retval==SYMSOLVER_WRONG_INERTIA &&
499                 augSysSolver_->NumberOfNegEVals() < numberOfEVals) {
500          Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
501                         "Number of negative eigenvalues too small!\n");
502          // If the number of negative eigenvalues is too small, then
503          // we first try to remedy this by asking for better quality
504          // solution (e.g. increasing pivot tolerance), and if that
505          // doesn't help, we assume that the system is singular
506          bool assume_singular = true;
507          if (!augsys_improved_) {
508            Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
509                           "Asking augmented system solver to improve quality of its solutions.\n");
510            augsys_improved_ = augSysSolver_->IncreaseQuality();
511            if (augsys_improved_) {
512              IpData().Append_info_string("q");
513              resolve_with_better_quality = true;
514              assume_singular = false;
515            }
516            else {
517              Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
518                             "Quality could not be improved\n");
519            }
520          }
521          if (assume_singular) {
522            bool pert_return =
523                               perturbHandler_->PerturbForSingularity(delta_x, delta_s,
524                                                                      delta_c, delta_d);
525            if (!pert_return) {
526              Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
527                             "PerturbForSingularity can't be done for assume singular.\n");
528              IpData().TimingStats().PDSystemSolverSolveOnce().End();
529              return false;
530            }
531            IpData().Append_info_string("a");
532          }
533        }
534        else if (retval==SYMSOLVER_WRONG_INERTIA ||
535                 retval==SYMSOLVER_SINGULAR) {
536          // Get new perturbation factors from the perturbation
537          // handlers for the case of wrong inertia
538          bool pert_return =
539                             perturbHandler_->PerturbForWrongInertia(delta_x, delta_s,
540                                                                     delta_c, delta_d);
541          if (!pert_return) {
542            Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
543                           "PerturbForWrongInertia can't be done for assume singular.\n");
544            IpData().TimingStats().PDSystemSolverSolveOnce().End();
545            return false;
546          }
547        }
548      } // while (retval!=SYMSOLVER_SUCCESS && !fail) {
549
550      // Some output
551      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
552                     "Number of trial factorizations performed: %d\n",
553                     count);
554      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
555                     "Perturbation parameters: delta_x=%e delta_s=%e\n                         delta_c=%e delta_d=%e\n",
556                     delta_x, delta_s, delta_c, delta_d);
557    }
558
559    // Compute the remaining sol Vectors
560    Px_L.SinvBlrmZMTdBr(-1., slack_x_L, *rhs.z_L(), z_L, *sol->x_NonConst(), *sol->z_L_NonConst());
561    Px_U.SinvBlrmZMTdBr(1., slack_x_U, *rhs.z_U(), z_U, *sol->x_NonConst(), *sol->z_U_NonConst());
562    Pd_L.SinvBlrmZMTdBr(-1., slack_s_L, *rhs.v_L(), v_L, *sol->s_NonConst(), *sol->v_L_NonConst());
563    Pd_U.SinvBlrmZMTdBr(1., slack_s_U, *rhs.v_U(), v_U, *sol->s_NonConst(), *sol->v_U_NonConst());
564
565    // Finally let's assemble the res result vectors
566    res.AddOneVector(alpha, *sol, beta);
567
568    IpData().TimingStats().PDSystemSolverSolveOnce().End();
569
570    return true;
571  }
572
573  void PDFullSpaceSolver::ComputeResiduals(
574    const SymMatrix& W,
575    const Matrix& J_c,
576    const Matrix& J_d,
577    const Matrix& Px_L,
578    const Matrix& Px_U,
579    const Matrix& Pd_L,
580    const Matrix& Pd_U,
581    const Vector& z_L,
582    const Vector& z_U,
583    const Vector& v_L,
584    const Vector& v_U,
585    const Vector& slack_x_L,
586    const Vector& slack_x_U,
587    const Vector& slack_s_L,
588    const Vector& slack_s_U,
589    const Vector& sigma_x,
590    const Vector& sigma_s,
591    Number alpha,
592    Number beta,
593    const IteratesVector& rhs,
594    const IteratesVector& res,
595    IteratesVector& resid)
596  {
597    DBG_START_METH("PDFullSpaceSolver::ComputeResiduals", dbg_verbosity);
598
599    DBG_PRINT_VECTOR(2, "res", res);
600    IpData().TimingStats().ComputeResiduals().Start();
601
602    // Get the current sizes of the perturbation factors
603    Number delta_x;
604    Number delta_s;
605    Number delta_c;
606    Number delta_d;
607    perturbHandler_->CurrentPerturbation(delta_x, delta_s, delta_c, delta_d);
608
609    SmartPtr<Vector> tmp;
610
611    // x
612    W.MultVector(1., *res.x(), 0., *resid.x_NonConst());
613    J_c.TransMultVector(1., *res.y_c(), 1., *resid.x_NonConst());
614    J_d.TransMultVector(1., *res.y_d(), 1., *resid.x_NonConst());
615    Px_L.MultVector(-1., *res.z_L(), 1., *resid.x_NonConst());
616    Px_U.MultVector(1., *res.z_U(), 1., *resid.x_NonConst());
617    resid.x_NonConst()->AddTwoVectors(delta_x, *res.x(), -1., *rhs.x(), 1.);
618
619    // s
620    Pd_U.MultVector(1., *res.v_U(), 0., *resid.s_NonConst());
621    Pd_L.MultVector(-1., *res.v_L(), 1., *resid.s_NonConst());
622    resid.s_NonConst()->AddTwoVectors(-1., *res.y_d(), -1., *rhs.s(), 1.);
623    if (delta_s!=0.) {
624      resid.s_NonConst()->Axpy(delta_s, *res.s());
625    }
626
627    // c
628    J_c.MultVector(1., *res.x(), 0., *resid.y_c_NonConst());
629    resid.y_c_NonConst()->AddTwoVectors(-delta_c, *res.y_c(), -1., *rhs.y_c(), 1.);
630
631    // d
632    J_d.MultVector(1., *res.x(), 0., *resid.y_d_NonConst());
633    resid.y_d_NonConst()->AddTwoVectors(-1., *res.s(), -1., *rhs.y_d(), 1.);
634    if (delta_d!=0.) {
635      resid.y_d_NonConst()->Axpy(-delta_d, *res.y_d());
636    }
637
638    // zL
639    resid.z_L_NonConst()->Copy(*res.z_L());
640    resid.z_L_NonConst()->ElementWiseMultiply(slack_x_L);
641    tmp = z_L.MakeNew();
642    Px_L.TransMultVector(1., *res.x(), 0., *tmp);
643    tmp->ElementWiseMultiply(z_L);
644    resid.z_L_NonConst()->AddTwoVectors(1., *tmp, -1., *rhs.z_L(), 1.);
645
646    // zU
647    resid.z_U_NonConst()->Copy(*res.z_U());
648    resid.z_U_NonConst()->ElementWiseMultiply(slack_x_U);
649    tmp = z_U.MakeNew();
650    Px_U.TransMultVector(1., *res.x(), 0., *tmp);
651    tmp->ElementWiseMultiply(z_U);
652    resid.z_U_NonConst()->AddTwoVectors(-1., *tmp, -1., *rhs.z_U(), 1.);
653
654    // vL
655    resid.v_L_NonConst()->Copy(*res.v_L());
656    resid.v_L_NonConst()->ElementWiseMultiply(slack_s_L);
657    tmp = v_L.MakeNew();
658    Pd_L.TransMultVector(1., *res.s(), 0., *tmp);
659    tmp->ElementWiseMultiply(v_L);
660    resid.v_L_NonConst()->AddTwoVectors(1., *tmp, -1., *rhs.v_L(), 1.);
661
662    // vU
663    resid.v_U_NonConst()->Copy(*res.v_U());
664    resid.v_U_NonConst()->ElementWiseMultiply(slack_s_U);
665    tmp = v_U.MakeNew();
666    Pd_U.TransMultVector(1., *res.s(), 0., *tmp);
667    tmp->ElementWiseMultiply(v_U);
668    resid.v_U_NonConst()->AddTwoVectors(-1., *tmp, -1., *rhs.v_U(), 1.);
669
670    DBG_PRINT_VECTOR(2, "resid", resid);
671
672    if (Jnlst().ProduceOutput(J_MOREVECTOR, J_LINEAR_ALGEBRA)) {
673      resid.Print(Jnlst(), J_MOREVECTOR, J_LINEAR_ALGEBRA, "resid");
674    }
675
676    if (Jnlst().ProduceOutput(J_MOREDETAILED, J_LINEAR_ALGEBRA)) {
677      Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
678                     "max-norm resid_x  %e\n", resid.x()->Amax());
679      Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
680                     "max-norm resid_s  %e\n", resid.s()->Amax());
681      Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
682                     "max-norm resid_c  %e\n", resid.y_c()->Amax());
683      Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
684                     "max-norm resid_d  %e\n", resid.y_d()->Amax());
685      Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
686                     "max-norm resid_zL %e\n", resid.z_L()->Amax());
687      Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
688                     "max-norm resid_zU %e\n", resid.z_U()->Amax());
689      Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
690                     "max-norm resid_vL %e\n", resid.v_L()->Amax());
691      Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
692                     "max-norm resid_vU %e\n", resid.v_U()->Amax());
693    }
694    IpData().TimingStats().ComputeResiduals().End();
695  }
696
697  Number PDFullSpaceSolver::ComputeResidualRatio(const IteratesVector& rhs,
698      const IteratesVector& res,
699      const IteratesVector& resid)
700  {
701    DBG_START_METH("PDFullSpaceSolver::ComputeResidualRatio", dbg_verbosity);
702
703    Number nrm_rhs = rhs.Amax();
704    Number nrm_res = res.Amax();
705    Number nrm_resid = resid.Amax();
706    Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
707                   "nrm_rhs = %8.2e nrm_sol = %8.2e nrm_resid = %8.2e\n",
708                   nrm_rhs, nrm_res, nrm_resid);
709
710    if (nrm_rhs+nrm_res == 0.) {
711      return nrm_resid;  // this should be zero
712    }
713    else {
714      // ToDo: determine how to include norm of matrix, and what
715      // safeguard to use against incredibly large solution vectors
716      Number max_cond = 1e6;
717      return nrm_resid/(Min(nrm_res, max_cond*nrm_rhs)+nrm_rhs);
718    }
719  }
720
721} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.