source: branches/dev/Algorithm/LinearSolvers/IpTSymLinearSolver.cpp @ 523

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

added TiminigStatistics? for measuring CPU time usage in different parts of
the algorithm. This might not yet compile on all platforms!
  • 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 523 2005-09-20 18:02:12Z 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.