source: branches/dev/Algorithm/LinearSolvers/IpTSymLinearSolver.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: 12.3 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: IpTSymLinearSolver.cpp 529 2005-09-29 21:12:38Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-03-17
8
9#include "IpTSymLinearSolver.hpp"
10#include "IpTripletHelper.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
22namespace Ipopt
23{
24#ifdef IP_DEBUG
25  static const Index dbg_verbosity = 0;
26#endif
27
28  TSymLinearSolver::TSymLinearSolver
29  (SmartPtr<SparseSymLinearSolverInterface> solver_interface,
30   SmartPtr<TSymScalingMethod> scaling_method)
31      :
32      SymLinearSolver(),
33      atag_(0),
34      dim_(0),
35      nonzeros_triplet_(0),
36      nonzeros_compressed_(0),
37      have_structure_(false),
38      initialized_(false),
39
40      solver_interface_(solver_interface),
41      scaling_method_(scaling_method),
42      scaling_factors_(NULL),
43      airn_(NULL),
44      ajcn_(NULL)
45  {
46    DBG_START_METH("TSymLinearSolver::TSymLinearSolver()",dbg_verbosity);
47    DBG_ASSERT(IsValid(solver_interface));
48  }
49
50  TSymLinearSolver::~TSymLinearSolver()
51  {
52    DBG_START_METH("TSymLinearSolver::~TSymLinearSolver()",
53                   dbg_verbosity);
54    delete [] airn_;
55    delete [] ajcn_;
56    delete [] scaling_factors_;
57  }
58
59  bool TSymLinearSolver::InitializeImpl(const OptionsList& options,
60                                        const std::string& prefix)
61  {
62    // This option is registered by OrigIpoptNLP
63    options.GetBoolValue("warm_start_same_structure",
64                         warm_start_same_structure_, prefix);
65
66    if (!solver_interface_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
67                                       options, prefix)) {
68      return false;
69    }
70
71    if (!warm_start_same_structure_) {
72      // Reset all private data
73      atag_=0;
74      dim_=0;
75      nonzeros_triplet_=0;
76      nonzeros_compressed_=0;
77      have_structure_=false;
78
79      matrix_format_ = solver_interface_->MatrixFormat();
80      switch (matrix_format_) {
81        case SparseSymLinearSolverInterface::CSR_Format_0_Offset:
82        triplet_to_csr_converter_ = new TripletToCSRConverter(0);
83        break;
84        case SparseSymLinearSolverInterface::CSR_Format_1_Offset:
85        triplet_to_csr_converter_ = new TripletToCSRConverter(1);
86        break;
87        case SparseSymLinearSolverInterface::Triplet_Format:
88        triplet_to_csr_converter_ = NULL;
89        break;
90        default:
91        DBG_ASSERT(false && "Invalid MatrixFormat returned from solver interface.");
92        return false;
93      }
94    }
95    else {
96      ASSERT_EXCEPTION(have_structure_, INVALID_WARMSTART,
97                       "TSymLinearSolver called with warm_start_same_structure, but the internal structures are not initialized.");
98    }
99
100    // reset the initialize flag to make sure that InitializeStructure
101    // is called for the linear solver
102    initialized_=false;
103
104    bool retval = true;
105    if (IsValid(scaling_method_)) {
106      IpData().TimingStats().LinearSystemScaling().Start();
107      retval = scaling_method_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
108                                           options, prefix);
109      IpData().TimingStats().LinearSystemScaling().End();
110    }
111    return retval;
112  }
113
114  ESymSolverStatus
115  TSymLinearSolver::MultiSolve(const SymMatrix& sym_A,
116                               std::vector<const Vector*>& rhsV,
117                               std::vector<Vector*>& solV,
118                               bool check_NegEVals,
119                               Index numberOfNegEVals)
120  {
121    DBG_START_METH("TSymLinearSolver::MultiSolve",dbg_verbosity);
122    DBG_ASSERT(!check_NegEVals || ProvidesInertia());
123
124    // Check if this object has ever seen a matrix If not,
125    // allocate memory of the matrix structure and copy the nonzeros
126    // structure (it is assumed that this will never change).
127    if (!initialized_) {
128      ESymSolverStatus retval = InitializeStructure(sym_A);
129      if (retval != SYMSOLVER_SUCCESS) {
130        return retval;
131      }
132    }
133
134    DBG_ASSERT(nonzeros_triplet_== TripletHelper::GetNumberEntries(sym_A));
135
136    // Check if the matrix has been changed
137    DBG_PRINT((1,"atag_=%d sym_A->GetTag()=%d\n",atag_,sym_A.GetTag()));
138    bool new_matrix = sym_A.HasChanged(atag_);
139    atag_ = sym_A.GetTag();
140
141    // If a new matrix is encountered, get the array for storing the
142    // entries from the linear solver interface, fill in the new
143    // values, compute the new scaling factors (if required), and
144    // scale the matrix
145    if (new_matrix) {
146      GiveMatrixToSolver(true, sym_A);
147    }
148
149    // Retrieve the right hand sides and scale if required
150    Index nrhs = (Index)rhsV.size();
151    double* rhs_vals = new double[dim_*nrhs];
152    for (Index irhs=0; irhs<nrhs; irhs++) {
153      TripletHelper::FillValuesFromVector(dim_, *rhsV[irhs],
154                                          &rhs_vals[irhs*(dim_)]);
155      if (IsValid(scaling_method_)) {
156        IpData().TimingStats().LinearSystemScaling().Start();
157        for (Index i=0; i<dim_; i++) {
158          rhs_vals[irhs*(dim_)+i] *= scaling_factors_[i];
159        }
160        IpData().TimingStats().LinearSystemScaling().End();
161      }
162    }
163
164    bool done = false;
165    // Call the linear solver through the interface to solve the
166    // system.  This is repeated, if the return values is S_CALL_AGAIN
167    // after the values have been restored (this might be necessary
168    // for MA27 if the size of the work space arrays was not large
169    // enough).
170    ESymSolverStatus retval;
171    while (!done) {
172      const Index* ia;
173      const Index* ja;
174      if (matrix_format_==SparseSymLinearSolverInterface::Triplet_Format) {
175        ia = airn_;
176        ja = ajcn_;
177      }
178      else {
179        ia = triplet_to_csr_converter_->IA();
180        ja = triplet_to_csr_converter_->JA();
181      }
182
183      retval = solver_interface_->MultiSolve(new_matrix, ia, ja,
184                                             nrhs, rhs_vals, check_NegEVals,
185                                             numberOfNegEVals);
186      if (retval==SYMSOLVER_CALL_AGAIN) {
187        DBG_PRINT((1, "Solver interface asks to be called again.\n"));
188        GiveMatrixToSolver(false, sym_A);
189      }
190      else {
191        done = true;
192      }
193    }
194
195    // If the solve was successful, unscale the solution (if required)
196    // and transfer the result into the Vectors
197    if (retval==SYMSOLVER_SUCCESS) {
198      for (Index irhs=0; irhs<nrhs; irhs++) {
199        if (IsValid(scaling_method_)) {
200          IpData().TimingStats().LinearSystemScaling().Start();
201          for (Index i=0; i<dim_; i++) {
202            rhs_vals[irhs*(dim_)+i] *= scaling_factors_[i];
203          }
204          IpData().TimingStats().LinearSystemScaling().End();
205        }
206        TripletHelper::PutValuesInVector(dim_, &rhs_vals[irhs*(dim_)],
207                                         *solV[irhs]);
208      }
209    }
210
211    delete[] rhs_vals;
212
213    return retval;
214  }
215
216  // Initialize the local copy of the positions of the nonzero
217  // elements
218  ESymSolverStatus
219  TSymLinearSolver::InitializeStructure(const SymMatrix& sym_A)
220  {
221    DBG_START_METH("TSymLinearSolver::InitializeStructure",
222                   dbg_verbosity);
223    DBG_ASSERT(!initialized_);
224
225    ESymSolverStatus retval;
226
227    // have_structure_ is already true if this is a warm start for a
228    // problem with identical structure
229    if (!have_structure_) {
230
231      dim_ = sym_A.Dim();
232      nonzeros_triplet_ = TripletHelper::GetNumberEntries(sym_A);
233
234      delete [] airn_;
235      delete [] ajcn_;
236      airn_ = new Index[nonzeros_triplet_];
237      ajcn_ = new Index[nonzeros_triplet_];
238
239      TripletHelper::FillRowCol(nonzeros_triplet_, sym_A, airn_, ajcn_);
240
241      // If the solver wants the compressed format, the converter has to
242      // be initialized
243      const Index *ia;
244      const Index *ja;
245      Index nonzeros;
246      if (matrix_format_ == SparseSymLinearSolverInterface::Triplet_Format) {
247        ia = airn_;
248        ja = ajcn_;
249        nonzeros = nonzeros_triplet_;
250      }
251      else {
252        nonzeros_compressed_ =
253          triplet_to_csr_converter_->InitializeConverter(dim_, nonzeros_triplet_,
254              airn_, ajcn_);
255        ia = triplet_to_csr_converter_->IA();
256        ja = triplet_to_csr_converter_->JA();
257        nonzeros = nonzeros_compressed_;
258      }
259
260      retval = solver_interface_->InitializeStructure(dim_, nonzeros, ia, ja);
261      if (retval != SYMSOLVER_SUCCESS) {
262        return retval;
263      }
264
265      // Get space for the scaling factors
266      delete [] scaling_factors_;
267      if (IsValid(scaling_method_)) {
268        IpData().TimingStats().LinearSystemScaling().Start();
269        scaling_factors_ = new double[dim_];
270        IpData().TimingStats().LinearSystemScaling().End();
271      }
272
273      have_structure_ = true;
274    }
275    else {
276      ASSERT_EXCEPTION(dim_==sym_A.Dim(), INVALID_WARMSTART,
277                       "TSymLinearSolver called with warm_start_same_structure, but the problem is solved for the first time.");
278      // This is a warm start for identical structure, so we don't need to
279      // recompute the nonzeros location arrays
280      const Index *ia;
281      const Index *ja;
282      Index nonzeros;
283      if (matrix_format_ == SparseSymLinearSolverInterface::Triplet_Format) {
284        ia = airn_;
285        ja = ajcn_;
286        nonzeros = nonzeros_triplet_;
287      }
288      else {
289        ia = triplet_to_csr_converter_->IA();
290        ja = triplet_to_csr_converter_->JA();
291        nonzeros = nonzeros_compressed_;
292      }
293      retval = solver_interface_->InitializeStructure(dim_, nonzeros, ia, ja);
294    }
295    initialized_=true;
296    return retval;
297  }
298
299  Index TSymLinearSolver::NumberOfNegEVals() const
300  {
301    DBG_START_METH("TSymLinearSolver::NumberOfNegEVals",dbg_verbosity);
302    return solver_interface_->NumberOfNegEVals();
303  }
304
305  bool TSymLinearSolver::IncreaseQuality()
306  {
307    DBG_START_METH("TSymLinearSolver::IncreaseQuality",dbg_verbosity);
308
309    return solver_interface_->IncreaseQuality();
310  }
311
312  bool TSymLinearSolver::ProvidesInertia() const
313  {
314    DBG_START_METH("TSymLinearSolver::ProvidesInertia",dbg_verbosity);
315
316    return solver_interface_->ProvidesInertia();
317  }
318
319  void TSymLinearSolver::GiveMatrixToSolver(bool new_matrix,
320      const SymMatrix& sym_A)
321  {
322    DBG_START_METH("TSymLinearSolver::GiveMatrixToSolver",dbg_verbosity);
323    DBG_PRINT((1,"new_matrix = %d\n",new_matrix));
324
325    double* pa = solver_interface_->GetValuesArrayPtr();
326    double* atriplet;
327
328    if (matrix_format_!=SparseSymLinearSolverInterface::Triplet_Format) {
329      atriplet = new double[nonzeros_triplet_];
330    }
331    else {
332      atriplet = pa;
333    }
334
335    //DBG_PRINT_MATRIX(3, "Aunscaled", sym_A);
336    TripletHelper::FillValues(nonzeros_triplet_, sym_A, atriplet);
337    if (DBG_VERBOSITY()>=3) {
338      for (Index i=0; i<nonzeros_triplet_; i++) {
339        DBG_PRINT((3, "KKTunscaled(%6d,%6d) = %24.16e\n", airn_[i], ajcn_[i], atriplet[i]));
340      }
341    }
342
343    if (IsValid(scaling_method_)) {
344      IpData().TimingStats().LinearSystemScaling().Start();
345      DBG_ASSERT(scaling_factors_);
346      if (new_matrix) {
347        // only compute scaling factors if the matrix has not been
348        // changed since the last call to this method
349        bool retval =
350          scaling_method_->ComputeSymTScalingFactors(dim_, nonzeros_triplet_,
351              airn_, ajcn_,
352              atriplet, scaling_factors_);
353        DBG_ASSERT(retval);
354        retval = false; // is added to make sure compiles doesn't
355        // complain if not in debug mode
356        if (Jnlst().ProduceOutput(J_MOREVECTOR, J_LINEAR_ALGEBRA)) {
357          for (Index i=0; i<dim_; i++) {
358            Jnlst().Printf(J_MOREVECTOR, J_LINEAR_ALGEBRA,
359                           "scaling factor[%6d] = %22.17e\n",
360                           i, scaling_factors_[i]);
361          }
362        }
363      }
364      for (Index i=0; i<nonzeros_triplet_; i++) {
365        atriplet[i] *=
366          scaling_factors_[airn_[i]-1] * scaling_factors_[ajcn_[i]-1];
367      }
368      if (DBG_VERBOSITY()>=3) {
369        for (Index i=0; i<nonzeros_triplet_; i++) {
370          DBG_PRINT((3, "KKTscaled(%6d,%6d) = %24.16e\n", airn_[i], ajcn_[i], atriplet[i]));
371        }
372      }
373      IpData().TimingStats().LinearSystemScaling().End();
374    }
375
376    if (matrix_format_!=SparseSymLinearSolverInterface::Triplet_Format) {
377      triplet_to_csr_converter_->ConvertValues(nonzeros_triplet_, atriplet,
378          nonzeros_compressed_, pa);
379      delete[] atriplet;
380    }
381
382  }
383
384} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.