source: branches/dev/Algorithm/IpStdAugSystemSolver.cpp @ 517

Last change on this file since 517 was 517, checked in by andreasw, 15 years ago
  • enabled warm start option for using same application and also for solving problems with identical structure
  • minor cosmetic changes
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.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: IpStdAugSystemSolver.cpp 517 2005-09-13 22:51:57Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8
9#include "IpStdAugSystemSolver.hpp"
10#include "IpDebug.hpp"
11
12#include "IpCompoundSymMatrix.hpp"
13#include "IpCompoundVector.hpp"
14#include "IpSumSymMatrix.hpp"
15#include "IpDiagMatrix.hpp"
16#include "IpIdentityMatrix.hpp"
17// ToDo: Remove below here - for debug only
18#include "IpTripletHelper.hpp"
19// ToDo: Remove above here
20
21namespace Ipopt
22{
23#ifdef IP_DEBUG
24  static const Index dbg_verbosity = 0;
25#endif
26
27  StdAugSystemSolver::StdAugSystemSolver(SymLinearSolver& linSolver)
28      :
29      AugSystemSolver(),
30      linsolver_(&linSolver),
31      augmented_system_space_(NULL),
32      sumsym_space_x_(NULL),
33      diag_space_x_(NULL),
34      diag_space_s_(NULL),
35      diag_space_c_(NULL),
36      ident_space_ds_(NULL),
37      diag_space_d_(NULL),
38      w_tag_(0),
39      d_x_tag_(0),
40      delta_x_(0.),
41      d_s_tag_(0),
42      delta_s_(0.),
43      j_c_tag_(0),
44      d_c_tag_(0),
45      delta_c_(0.),
46      j_d_tag_(0),
47      d_d_tag_(0),
48      delta_d_(0.),
49      old_w_(NULL)
50  {
51    DBG_START_METH("StdAugSystemSolver::StdAugSystemSolver()",dbg_verbosity);
52    DBG_ASSERT(IsValid(linsolver_));
53  }
54
55  StdAugSystemSolver::~StdAugSystemSolver()
56  {
57    DBG_START_METH("StdAugSystemSolver::~StdAugSystemSolver()",dbg_verbosity);
58  }
59
60
61  bool StdAugSystemSolver::InitializeImpl(const OptionsList& options,
62                                          const std::string& prefix)
63  {
64    // This option is registered by OrigIpoptNLP
65    options.GetBoolValue("warm_start_same_structure",
66                         warm_start_same_structure_, prefix);
67
68    if (!warm_start_same_structure_) {
69      augsys_tag_ = 0;
70      augmented_system_ = NULL;
71    }
72    else {
73      ASSERT_EXCEPTION(IsValid(augmented_system_), INVALID_WARMSTART,
74                       "StdAugSystemSolver called with warm_start_same_structure, but augmented system is not initialized.");
75    }
76
77    return linsolver_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
78                                  options, prefix);
79  }
80
81
82  ESymSolverStatus StdAugSystemSolver::Solve(const SymMatrix* W,
83      const Vector* D_x,
84      double delta_x,
85      const Vector* D_s,
86      double delta_s,
87      const Matrix* J_c,
88      const Vector* D_c,
89      double delta_c,
90      const Matrix* J_d,
91      const Vector* D_d,
92      double delta_d,
93      const Vector& rhs_x,
94      const Vector& rhs_s,
95      const Vector& rhs_c,
96      const Vector& rhs_d,
97      Vector& sol_x,
98      Vector& sol_s,
99      Vector& sol_c,
100      Vector& sol_d,
101      bool check_NegEVals,
102      Index numberOfNegEVals)
103  {
104    DBG_START_METH("StdAugSystemSolver::Solve",dbg_verbosity);
105    DBG_ASSERT(J_c && J_d && "Currently, you MUST specify J_c and J_d in the augmented system");
106
107    // Create the compound matrix of the augmented system if it has not
108    // yet been created - It is assumed that the structure will not change
109    // after this call
110    bool debug_first_time_through = false;
111    if (!IsValid(augmented_system_)) {
112      // pass in the information to form the structure of the augmented system
113      // rhs_? are passed in to provide a prototype vector
114      // for D_? (since these may be NULL)
115      DBG_ASSERT(W && J_c && J_d); // W must exist during the first call to setup the structure!
116      CreateAugmentedSpace(*W, *J_c, *J_d, rhs_x, rhs_s, rhs_c, rhs_d);
117      CreateAugmentedSystem(W, D_x, delta_x, D_s, delta_s,
118                            *J_c, D_c, delta_c, *J_d, D_d, delta_d,
119                            rhs_x, rhs_s, rhs_c, rhs_d);
120      debug_first_time_through = true;
121    }
122
123
124    // Check if anything that was just passed in is different from what is currently
125    // in the compound matrix of the augmented system. If anything is different, then
126    // update the augmented system
127    if ( AugmentedSystemRequiresChange(W, D_x, delta_x, D_s, delta_s, *J_c, D_c, delta_c, *J_d, D_d, delta_d) ) {
128      DBG_ASSERT(!debug_first_time_through);
129      CreateAugmentedSystem(W, D_x, delta_x, D_s, delta_s,
130                            *J_c, D_c, delta_c, *J_d, D_d, delta_d,
131                            rhs_x, rhs_s, rhs_c, rhs_d);
132    }
133
134    // Sanity checks
135    DBG_ASSERT(rhs_x.Dim() == sol_x.Dim());
136    DBG_ASSERT(rhs_s.Dim() == sol_s.Dim());
137    DBG_ASSERT(rhs_c.Dim() == sol_c.Dim());
138    DBG_ASSERT(rhs_d.Dim() == sol_d.Dim());
139
140    // Now construct the overall right hand side vector that will be passed
141    // to the linear solver
142    SmartPtr<CompoundVector> augmented_rhs = augmented_vector_space_->MakeNewCompoundVector();
143    augmented_rhs->SetComp(0, rhs_x);
144    augmented_rhs->SetComp(1, rhs_s);
145    augmented_rhs->SetComp(2, rhs_c);
146    augmented_rhs->SetComp(3, rhs_d);
147
148    augmented_system_->Print(Jnlst(), J_MATRIX, J_LINEAR_ALGEBRA, "KKT");
149    if (Jnlst().ProduceOutput(J_MOREMATRIX, J_LINEAR_ALGEBRA)) {
150      // ToDo: remove below here - for debug only
151      Index dbg_nz = TripletHelper::GetNumberEntries(*augmented_system_);
152      Index* dbg_iRows = new Index[dbg_nz];
153      Index* dbg_jCols = new Index[dbg_nz];
154      Number* dbg_values = new Number[dbg_nz];
155      TripletHelper::FillRowCol(dbg_nz, *augmented_system_, dbg_iRows, dbg_jCols);
156      TripletHelper::FillValues(dbg_nz, *augmented_system_, dbg_values);
157      Jnlst().Printf(J_MOREMATRIX, J_LINEAR_ALGEBRA, "******* KKT SYSTEM *******\n");
158      for (Index dbg_i=0; dbg_i<dbg_nz; dbg_i++) {
159        Jnlst().Printf(J_MOREMATRIX, J_LINEAR_ALGEBRA, "(%d) KKT[%d][%d] = %23.15e\n", dbg_i, dbg_iRows[dbg_i], dbg_jCols[dbg_i], dbg_values[dbg_i]);
160      }
161      delete [] dbg_iRows;
162      dbg_iRows = NULL;
163      delete [] dbg_jCols;
164      dbg_jCols = NULL;
165      delete [] dbg_values;
166      dbg_values = NULL;
167      // ToDo: remove above here
168    }
169    augmented_rhs->Print(Jnlst(), J_MOREVECTOR, J_LINEAR_ALGEBRA, "RHS");
170
171    // Call the linear solver
172    SmartPtr<CompoundVector> augmented_sol = augmented_vector_space_->MakeNewCompoundVector();
173    augmented_sol->SetCompNonConst(0, sol_x);
174    augmented_sol->SetCompNonConst(1, sol_s);
175    augmented_sol->SetCompNonConst(2, sol_c);
176    augmented_sol->SetCompNonConst(3, sol_d);
177    ESymSolverStatus retval;
178    retval = linsolver_->Solve(*augmented_system_, *augmented_rhs, *augmented_sol,
179                               check_NegEVals, numberOfNegEVals);
180    if (retval==SYMSOLVER_SUCCESS) {
181      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "Factorization successful.\n");
182      augmented_sol->Print(Jnlst(), J_MOREVECTOR, J_LINEAR_ALGEBRA, "SOL");
183    }
184    else if (retval==SYMSOLVER_FATAL_ERROR) {
185      THROW_EXCEPTION(FATAL_ERROR_IN_LINEAR_SOLVER,"A fatal error occured in the linear solver.");
186    }
187    else {
188      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "Factorization failed with retval = %d\n", retval);
189    }
190
191    return retval;
192  }
193
194  void StdAugSystemSolver::CreateAugmentedSpace(const SymMatrix& W,
195      const Matrix& J_c,
196      const Matrix& J_d,
197      const Vector& proto_x,
198      const Vector& proto_s,
199      const Vector& proto_c,
200      const Vector& proto_d)
201  {
202    DBG_ASSERT(!IsValid(augmented_system_));
203
204    old_w_ = &W;
205
206    //===
207    // Setup the augmented system matrix (described in IpAugSystemSolver.hpp")
208    //===
209
210    // created the compound symmetric matrix space
211    Index n_x = J_c.NCols();
212    Index n_s = J_d.NRows();
213    Index n_c = J_c.NRows();
214    Index n_d = n_s;
215
216    Index total_nRows = n_x + n_s + n_c + n_d;
217    augmented_system_space_ = new CompoundSymMatrixSpace(4, total_nRows);
218    augmented_system_space_->SetBlockDim(0, n_x);
219    augmented_system_space_->SetBlockDim(1, n_s);
220    augmented_system_space_->SetBlockDim(2, n_c);
221    augmented_system_space_->SetBlockDim(3, n_d);
222
223    // (1,1) block
224    // create the spaces and sum matrix for the upper left corner (W + D_x delta_x*I)
225    // of the hessian part for the 1,1 block
226    sumsym_space_x_ = new SumSymMatrixSpace(n_x, 2);
227    augmented_system_space_->SetCompSpace(0,0, *sumsym_space_x_);
228
229    diag_space_x_ = new DiagMatrixSpace(n_x);
230
231    // (2,2) block
232    // create the spaces and diag matrix for the lower right corner (D_s + delta_s*I)
233    // of the hessian part, the 2,2 block
234    diag_space_s_ = new DiagMatrixSpace(n_s);
235    augmented_system_space_->SetCompSpace(1,1, *diag_space_s_);
236
237    // (3,1) block
238    augmented_system_space_->SetCompSpace(2,0, *J_c.OwnerSpace());
239
240    // (3,3) block
241    // create the matrix space and matrix for the 3,3 block
242    diag_space_c_ = new DiagMatrixSpace(n_c);
243    augmented_system_space_->SetCompSpace(2,2, *diag_space_c_);
244
245    // (4,1) block
246    augmented_system_space_->SetCompSpace(3,0, *J_d.OwnerSpace());
247
248    // (4,2) block
249    // create the identity matrix space and matrix for the 4,2 block
250    ident_space_ds_ = new IdentityMatrixSpace(n_s);
251    augmented_system_space_->SetCompSpace(3,1, *ident_space_ds_);
252
253    // (4,4) block
254    // create the sum matrix space and matrix for the 4,4 block
255    diag_space_d_ = new DiagMatrixSpace(n_d);
256    augmented_system_space_->SetCompSpace(3,3, *diag_space_d_);
257
258    // Create the space for the vectors
259    augmented_vector_space_ = new CompoundVectorSpace(4, n_x + n_s + n_c + n_d);
260    augmented_vector_space_->SetCompSpace(0, *proto_x.OwnerSpace());
261    augmented_vector_space_->SetCompSpace(1, *proto_s.OwnerSpace());
262    augmented_vector_space_->SetCompSpace(2, *proto_c.OwnerSpace());
263    augmented_vector_space_->SetCompSpace(3, *proto_d.OwnerSpace());
264
265  }
266
267  void StdAugSystemSolver::CreateAugmentedSystem(const SymMatrix* W,
268      const Vector* D_x,
269      double delta_x,
270      const Vector* D_s,
271      double delta_s,
272      const Matrix& J_c,
273      const Vector* D_c,
274      double delta_c,
275      const Matrix& J_d,
276      const Vector* D_d,
277      double delta_d,
278      const Vector& proto_x,
279      const Vector& proto_s,
280      const Vector& proto_c,
281      const Vector& proto_d)
282  {
283    augmented_system_ = augmented_system_space_->MakeNewCompoundSymMatrix();
284
285    // (1,1) block
286    SmartPtr<SumSymMatrix> sumsym_x = sumsym_space_x_->MakeNewSumSymMatrix();
287
288    if (W) {
289      sumsym_x->SetTerm(0, 1.0, *W);
290      old_w_ = W;
291      w_tag_ = W->GetTag();
292    }
293    else {
294      sumsym_x->SetTerm(0, 0.0, *old_w_);
295      w_tag_ = 0;
296    }
297
298    SmartPtr<DiagMatrix> diag_x = diag_space_x_->MakeNewDiagMatrix();
299    if (D_x) {
300      if (delta_x==0.) {
301        diag_x->SetDiag(*D_x);
302      }
303      else {
304        SmartPtr<Vector> tmp = D_x->MakeNewCopy();
305        tmp->AddScalar(delta_x);
306        diag_x->SetDiag(*tmp);
307      }
308      d_x_tag_ = D_x->GetTag();
309    }
310    else {
311      SmartPtr<Vector> tmp = proto_x.MakeNew();
312      tmp->Set(delta_x);
313      diag_x->SetDiag(*tmp);
314      d_x_tag_ = 0;
315    }
316    sumsym_x->SetTerm(1, 1.0, *diag_x);
317    delta_x_ = delta_x;
318
319    augmented_system_->SetComp(0,0, *sumsym_x);
320
321    // (2,2) block
322    SmartPtr<DiagMatrix> diag_s = diag_space_s_->MakeNewDiagMatrix();
323    if (D_s) {
324      if (delta_s==0.) {
325        diag_s->SetDiag(*D_s);
326      }
327      else {
328        SmartPtr<Vector> tmp = D_s->MakeNewCopy();
329        tmp->AddScalar(delta_s);
330        diag_s->SetDiag(*tmp);
331      }
332      d_s_tag_ = D_s->GetTag();
333    }
334    else {
335      SmartPtr<Vector> tmp = proto_s.MakeNew();
336      tmp->Set(delta_s);
337      diag_s->SetDiag(*tmp);
338      d_s_tag_ = 0;
339    }
340    delta_s_ = delta_s;
341
342    augmented_system_->SetComp(1, 1, *diag_s);
343
344    // (3,1) block
345    augmented_system_->SetComp(2,0, J_c);
346    j_c_tag_ = J_c.GetTag();
347
348    // (3,3) block
349    SmartPtr<DiagMatrix> diag_c = diag_space_c_->MakeNewDiagMatrix();
350    if (D_c) {
351      if (delta_c==0.) {
352        diag_c->SetDiag(*D_c);
353      }
354      else {
355        SmartPtr<Vector> tmp = D_c->MakeNewCopy();
356        tmp->AddScalar(-delta_c);
357        diag_c->SetDiag(*tmp);
358      }
359      d_c_tag_ = D_c->GetTag();
360    }
361    else {
362      SmartPtr<Vector> tmp = proto_c.MakeNew();
363      tmp->Set(-delta_c);
364      diag_c->SetDiag(*tmp);
365      d_c_tag_ = 0;
366    }
367    delta_c_ = delta_c;
368
369    augmented_system_->SetComp(2,2, *diag_c);
370
371    // (4,1) block
372    augmented_system_->SetComp(3,0, J_d);
373    j_d_tag_ = J_d.GetTag();
374
375    // (4,2) block
376    SmartPtr<IdentityMatrix> ident_ds = ident_space_ds_->MakeNewIdentityMatrix();
377    ident_ds->SetFactor(-1.0);
378    augmented_system_->SetComp(3,1, *ident_ds);
379
380    // (4,4) block
381    SmartPtr<DiagMatrix> diag_d = diag_space_d_->MakeNewDiagMatrix();
382    if (D_d) {
383      if (delta_d==0.) {
384        diag_d->SetDiag(*D_d);
385      }
386      else {
387        SmartPtr<Vector> tmp = D_d->MakeNewCopy();
388        tmp->AddScalar(-delta_d);
389        diag_d->SetDiag(*tmp);
390      }
391      d_d_tag_ = D_d->GetTag();
392    }
393    else {
394      SmartPtr<Vector> tmp = proto_d.MakeNew();
395      tmp->Set(-delta_d);
396      diag_d->SetDiag(*tmp);
397      d_d_tag_ = 0;
398    }
399    delta_d_ = delta_d;
400
401    augmented_system_->SetComp(3,3, *diag_d);
402
403    augsys_tag_ = augmented_system_->GetTag();
404  }
405
406
407  bool StdAugSystemSolver::AugmentedSystemRequiresChange(const SymMatrix* W,
408      const Vector* D_x,
409      double delta_x,
410      const Vector* D_s,
411      double delta_s,
412      const Matrix& J_c,
413      const Vector* D_c,
414      double delta_c,
415      const Matrix& J_d,
416      const Vector* D_d,
417      double delta_d)
418  {
419    DBG_START_METH("StdAugSystemSolver::AugmentedSystemRequiresChange",dbg_verbosity);
420    DBG_ASSERT(augsys_tag_ == augmented_system_->GetTag() && "Someone has changed the augmented system outside of the AugSystemSolver. This should NOT happen.");
421
422#ifdef IP_DEBUG
423
424    bool Wtest = (W && W->GetTag() != w_tag_);
425    bool iWtest = (!W && w_tag_ != 0);
426    bool D_xtest = (D_x && D_x->GetTag() != d_x_tag_);
427    bool iD_xtest = (!D_x && d_x_tag_ != 0);
428    bool delta_xtest = (delta_x != delta_x_);
429    bool D_stest = (D_s && D_s->GetTag() != d_s_tag_);
430    bool iD_stest = (!D_s && d_s_tag_ != 0);
431    bool delta_stest = (delta_s != delta_s_);
432    bool J_ctest = (J_c.GetTag() != j_c_tag_);
433    bool D_ctest = (D_c && D_c->GetTag() != d_c_tag_);
434    bool iD_ctest = (!D_c && d_c_tag_ != 0);
435    bool delta_ctest = (delta_c != delta_c_);
436    bool J_dtest = (J_d.GetTag() != j_d_tag_);
437    bool D_dtest = (D_d && D_d->GetTag() != d_d_tag_);
438    bool iD_dtest = (!D_d && d_d_tag_ != 0);
439    bool delta_dtest = (delta_d != delta_d_);
440#endif
441
442    DBG_PRINT((2,"Wtest = %d\n", Wtest));
443    DBG_PRINT((2,"iWtest = %d\n", iWtest));
444    DBG_PRINT((2,"D_xtest = %d\n", D_xtest));
445    DBG_PRINT((2,"iD_xtest = %d\n", iD_xtest));
446    DBG_PRINT((2,"delta_xtest = %d\n", delta_xtest));
447    DBG_PRINT((2,"D_stest = %d\n", D_stest));
448    DBG_PRINT((2,"iD_stest = %d\n", iD_stest));
449    DBG_PRINT((2,"delta_stest = %d\n", delta_stest));
450    DBG_PRINT((2,"J_ctest = %d\n", J_ctest));
451    DBG_PRINT((2,"D_ctest = %d\n", D_ctest));
452    DBG_PRINT((2,"iD_ctest = %d\n", iD_ctest));
453    DBG_PRINT((2,"delta_ctest = %d\n", delta_ctest));
454    DBG_PRINT((2,"J_dtest = %d\n", J_dtest));
455    DBG_PRINT((2,"D_dtest = %d\n", D_dtest));
456    DBG_PRINT((2,"iD_dtest = %d\n", iD_dtest));
457    DBG_PRINT((2,"delta_dtest = %d\n", delta_dtest));
458
459    if ( (W && W->GetTag() != w_tag_)
460         || (!W && w_tag_ != 0)
461         || (D_x && D_x->GetTag() != d_x_tag_)
462         || (!D_x && d_x_tag_ != 0)
463         || (delta_x != delta_x_)
464         || (D_s && D_s->GetTag() != d_s_tag_)
465         || (!D_s && d_s_tag_ != 0)
466         || (delta_s != delta_s_)
467         || (J_c.GetTag() != j_c_tag_)
468         || (D_c && D_c->GetTag() != d_c_tag_)
469         || (!D_c && d_c_tag_ != 0)
470         || (delta_c != delta_c_)
471         || (J_d.GetTag() != j_d_tag_)
472         || (D_d && D_d->GetTag() != d_d_tag_)
473         || (!D_d && d_d_tag_ != 0)
474         || (delta_d != delta_d_) ) {
475      return true;
476    }
477
478    return false;
479  }
480
481  Index StdAugSystemSolver::NumberOfNegEVals() const
482  {
483    DBG_ASSERT(IsValid(augmented_system_));
484    return linsolver_->NumberOfNegEVals();
485  }
486
487  bool StdAugSystemSolver::ProvidesInertia() const
488  {
489    return linsolver_->ProvidesInertia();
490  }
491
492  bool StdAugSystemSolver::IncreaseQuality()
493  {
494    return linsolver_->IncreaseQuality();
495  }
496
497} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.