source: branches/dev/LinAlg/IpVector.hpp @ 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: 20.2 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: IpVector.hpp 529 2005-09-29 21:12:38Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8
9#ifndef __IPVECTOR_HPP__
10#define __IPVECTOR_HPP__
11
12#include "IpTypes.hpp"
13#include "IpTaggedObject.hpp"
14#include "IpCachedResults.hpp"
15#include "IpSmartPtr.hpp"
16#include "IpJournalist.hpp"
17
18#include <vector>
19
20#ifdef HAVE_CMATH
21# include <cmath>
22#else
23# ifdef HAVE_MATH_H
24#  include <math.h>
25# else
26#  error "don't have header file for math"
27# endif
28#endif
29
30namespace Ipopt
31{
32  /* forward declarations */
33  class VectorSpace;
34
35  /** Vector Base Class.
36   * This is the base class for all derived vector types.  Those vectors
37   * are meant to store entities like iterates, Lagrangian multipliers,
38   * constraint values etc.  The implementation of a vector type depends
39   * on the computational environment (e.g. just a double array on a shared
40   * memory machine, or distributed double arrays for a distributed
41   * memory machine.)
42   *
43   * Deriving from Vector: This class inherits from tagged object to
44   * implement an advanced caching scheme. Because of this, the
45   * TaggedObject method ObjectChanged() must be called each time the
46   * Vector changes. If you overload the XXXX_Impl protected methods,
47   * this taken care of (along with caching if possible) for you. If
48   * you have additional methods in your derived class that change the
49   * underlying data (vector values), you MUST remember to call
50   * ObjectChanged() AFTER making the change!
51   */
52  class Vector : public TaggedObject
53  {
54  public:
55    /** @name Constructor/Destructor */
56    //@{
57    /** Constructor.  It has to be given a pointer to the
58     *  corresponding VectorSpace.
59     */
60    Vector(const VectorSpace* owner_space);
61
62    /** Destructor */
63    virtual ~Vector()
64    {}
65    //@}
66
67    /** Create new Vector of the same type with uninitialized data */
68    Vector* MakeNew() const;
69
70    /** Create new Vector of the same type and copy the data over */
71    Vector* MakeNewCopy() const;
72
73    /**@name Standard BLAS-1 Operations
74     *  (derived classes do NOT overload these
75     *  methods, instead, overload the
76     *  protected versions of these methods). */
77    //@{
78    /** Copy the data of the vector x into this vector (DCOPY). */
79    void Copy(const Vector& x);
80
81    /** Scales the vector by scalar alpha (DSCAL) */
82    void Scal(Number alpha);
83
84    /** Add the multiple alpha of vector x to this vector (DAXPY) */
85    void Axpy(Number alpha, const Vector &x);
86
87    /** Computes inner product of vector x with this (DDOT) */
88    Number Dot(const Vector &x) const;
89
90    /** Computes the 2-norm of this vector (DNRM2) */
91    Number Nrm2() const;
92
93    /** Computes the 1-norm of this vector (DASUM) */
94    Number Asum() const;
95
96    /** Computes the max-norm of this vector (based on IDAMAX) */
97    Number Amax() const;
98    //@}
99
100    /** @name Additional (Non-BLAS) Vector Methods
101     *  (derived classes do NOT overload these
102     *  methods, instead, overload the
103     *  protected versions of these methods). */
104    //@{
105    /** Set each element in the vector to the scalar alpha. */
106    void Set(Number alpha);
107
108    /** Element-wise division  \f$y_i \gets y_i/x_i\f$*/
109    void ElementWiseDivide(const Vector& x);
110
111    /** Element-wise multiplication \f$y_i \gets y_i*x_i\f$ */
112    void ElementWiseMultiply(const Vector& x);
113
114    /** Element-wise max against entries in x */
115    void ElementWiseMax(const Vector& x);
116
117    /** Element-wise min against entries in x */
118    void ElementWiseMin(const Vector& x);
119
120    /** Reciprocates the entries in the vector */
121    void ElementWiseReciprocal();
122
123    /** Absolute values of the entries in the vector */
124    void ElementWiseAbs();
125
126    /** Element-wise square root of the entries in the vector */
127    void ElementWiseSqrt();
128
129    /** Replaces the vector values with their sgn values
130    ( -1 if x_i < 0, 0 if x_i == 0, and 1 if x_i > 0)
131    */
132    void ElementWiseSgn();
133
134    /** Add scalar to every vector component */
135    void AddScalar(Number scalar);
136
137    /** Returns the maximum value in the vector */
138    Number Max() const;
139
140    /** Returns the minimum value in the vector */
141    Number Min() const;
142
143    /** Returns the sum of the vector entries */
144    Number Sum() const;
145
146    /** Returns the sum of the logs of each vector entry */
147    Number SumLogs() const;
148    //@}
149
150    /** @name Methods for specialized operations.  A prototype
151     *  implementation is provided, but for efficient implementation
152     *  those should be specially implemented.
153     */
154    //@{
155    /** Add one vector, y = a * v1 + c * y.  This is automatically
156     *  reduced to call AddTwoVectors.  */
157    void AddOneVector(Number a, const Vector& v1, Number c);
158
159    /** Add two vectors, y = a * v1 + b * v2 + c * y.  Here, this
160     *  vector is y */
161    void AddTwoVectors(Number a, const Vector& v1,
162                       Number b, const Vector& v2, Number c);
163    /** Fraction to the boundary parameter.  Computes \f$\alpha =
164     *  \max\{\bar\alpha\in(0,1] : x + \bar\alpha \Delta \geq (1-\tau)x\}\f$
165     */
166    Number FracToBound(const Vector& delta, Number tau) const;
167    /** Add the quotient of two vectors, y = a * z/s + c * y. */
168    void AddVectorQuotient(Number a, const Vector& z, const Vector& s,
169                           Number c);
170    //@}
171
172    /** @name Accessor methods */
173    //@{
174    /** Dimension of the Vector */
175    Index Dim() const;
176
177    /** Return the owner VectorSpace*/
178    SmartPtr<const VectorSpace> OwnerSpace() const;
179    //@}
180
181    /** @name Output methods
182     *  (derived classes do NOT overload these
183     *  methods, instead, overload the
184     *  protected versions of these methods). */
185    //@{
186    /** Print the entire vector */
187    void Print(SmartPtr<const Journalist> jnlst,
188               EJournalLevel level,
189               EJournalCategory category,
190               const std::string& name,
191               Index indent=0,
192               const std::string& prefix="") const;
193    void Print(const Journalist& jnlst,
194               EJournalLevel level,
195               EJournalCategory category,
196               const std::string& name,
197               Index indent=0,
198               const std::string& prefix="") const;
199    //@}
200
201  protected:
202    /** @name implementation methods (derived classes MUST
203     *  overload these pure virtual protected methods.)
204     */
205    //@{
206    /** Copy the data of the vector x into this vector (DCOPY). */
207    virtual void CopyImpl(const Vector& x)=0;
208
209    /** Scales the vector by scalar alpha (DSCAL) */
210    virtual void ScalImpl(Number alpha)=0;
211
212    /** Add the multiple alpha of vector x to this vector (DAXPY) */
213    virtual void AxpyImpl(Number alpha, const Vector &x)=0;
214
215    /** Computes inner product of vector x with this (DDOT) */
216    virtual Number DotImpl(const Vector &x) const =0;
217
218    /** Computes the 2-norm of this vector (DNRM2) */
219    virtual Number Nrm2Impl() const =0;
220
221    /** Computes the 1-norm of this vector (DASUM) */
222    virtual Number AsumImpl() const =0;
223
224    /** Computes the max-norm of this vector (based on IDAMAX) */
225    virtual Number AmaxImpl() const =0;
226
227    /** Set each element in the vector to the scalar alpha. */
228    virtual void SetImpl(Number alpha)=0;
229
230    /** Element-wise division  \f$y_i \gets y_i/x_i\f$*/
231    virtual void ElementWiseDivideImpl(const Vector& x)=0;
232
233    /** Element-wise multiplication \f$y_i \gets y_i*x_i\f$ */
234    virtual void ElementWiseMultiplyImpl(const Vector& x)=0;
235
236    /** Element-wise max against entries in x */
237    virtual void ElementWiseMaxImpl(const Vector& x)=0;
238
239    /** Element-wise min against entries in x */
240    virtual void ElementWiseMinImpl(const Vector& x)=0;
241
242    /** Reciprocates the elements of the vector */
243    virtual void ElementWiseReciprocalImpl()=0;
244
245    /** Take elementwise absolute values of the elements of the vector */
246    virtual void ElementWiseAbsImpl()=0;
247
248    /** Take elementwise square-root of the elements of the vector */
249    virtual void ElementWiseSqrtImpl()=0;
250
251    /** Replaces entries with sgn of the entry */
252    virtual void ElementWiseSgnImpl()=0;
253
254    /** Add scalar to every component of vector */
255    virtual void AddScalarImpl(Number scalar)=0;
256
257    /** Max value in the vector */
258    virtual Number MaxImpl() const=0;
259
260    /** Min number in the vector */
261    virtual Number MinImpl() const=0;
262
263    /** Sum of entries in the vector */
264    virtual Number SumImpl() const=0;
265
266    /** Sum of logs of entries in the vector */
267    virtual Number SumLogsImpl() const=0;
268
269    /** Add two vectors (a * v1 + b * v2).  Result is stored in this
270    vector. */
271    virtual void AddTwoVectorsImpl(Number a, const Vector& v1,
272                                   Number b, const Vector& v2, Number c);
273
274    /** Fraction to boundary parameter. */
275    virtual Number FracToBoundImpl(const Vector& delta, Number tau) const;
276
277    /** Add the quotient of two vectors */
278    virtual void AddVectorQuotientImpl(Number a, const Vector& z,
279                                       const Vector& s, Number c);
280
281    /** Print the entire vector */
282    virtual void PrintImpl(const Journalist& jnlst,
283                           EJournalLevel level,
284                           EJournalCategory category,
285                           const std::string& name,
286                           Index indent,
287                           const std::string& prefix) const =0;
288    //@}
289
290  private:
291    /**@name Default Compiler Generated Methods
292     * (Hidden to avoid implicit creation/calling).
293     * These methods are not implemented and
294     * we do not want the compiler to implement
295     * them for us, so we declare them private
296     * and do not define them. This ensures that
297     * they will not be implicitly created/called. */
298    //@{
299    /** Default constructor */
300    Vector();
301
302    /** Copy constructor */
303    Vector(const Vector&);
304
305    /** Overloaded Equals Operator */
306    Vector& operator=(const Vector&);
307    //@}
308
309    /** Vector Space */
310    const SmartPtr<const VectorSpace> owner_space_;
311
312    /**@name CachedResults data members */
313    //@{
314    /** Cache for dot products */
315    mutable CachedResults<Number> dot_cache_;
316
317    mutable TaggedObject::Tag nrm2_cache_tag_;
318    mutable Number cached_nrm2_;
319
320    mutable TaggedObject::Tag asum_cache_tag_;
321    mutable Number cached_asum_;
322
323    mutable TaggedObject::Tag amax_cache_tag_;
324    mutable Number cached_amax_;
325
326    mutable TaggedObject::Tag max_cache_tag_;
327    mutable Number cached_max_;
328
329    mutable TaggedObject::Tag min_cache_tag_;
330    mutable Number cached_min_;
331
332    mutable TaggedObject::Tag sum_cache_tag_;
333    mutable Number cached_sum_;
334
335    mutable TaggedObject::Tag sumlogs_cache_tag_;
336    mutable Number cached_sumlogs_;
337
338    //     AW: I removed this cache since it gets in the way for the
339    //         quality function search
340    //     /** Cache for FracToBound */
341    //     mutable CachedResults<Number> frac_to_bound_cache_;
342    //@}
343
344  };
345
346  /** VectorSpace base class, corresponding to the Vector base class.
347   *  For each Vector implementation, a corresponding VectorSpace has
348   *  to be implemented.  A VectorSpace is able to create new Vectors
349   *  of a specific type.  The VectorSpace should also store
350   *  information that is common to all Vectors of that type.  For
351   *  example, the dimension of a Vector is stored in the VectorSpace
352   *  base class.
353   */
354  class VectorSpace : public ReferencedObject
355  {
356  public:
357    /** @name Constructors/Destructors */
358    //@{
359    /** Constructor, given the dimension of all vectors generated by
360     *  this VectorSpace.
361     */
362    VectorSpace(Index dim);
363
364    /** Destructor */
365    virtual ~VectorSpace()
366    {}
367    //@}
368
369    /** Pure virtual method for creating a new Vector of the
370     *  corresponding type.
371     */
372    virtual Vector* MakeNew() const=0;
373
374    /** Accessor function for the dimension of the vectors of this type.*/
375    Index Dim() const
376    {
377      return dim_;
378    }
379
380  private:
381    /**@name Default Compiler Generated Methods
382     * (Hidden to avoid implicit creation/calling).
383     * These methods are not implemented and
384     * we do not want the compiler to implement
385     * them for us, so we declare them private
386     * and do not define them. This ensures that
387     * they will not be implicitly created/called. */
388    //@{
389    /** default constructor */
390    VectorSpace();
391
392    /** Copy constructor */
393    VectorSpace(const VectorSpace&);
394
395    /** Overloaded Equals Operator */
396    VectorSpace& operator=(const VectorSpace&);
397    //@}
398
399    /** Dimension of the vectors in this vector space. */
400    const Index dim_;
401  };
402
403  /* inline methods */
404  inline
405  Vector::Vector(const VectorSpace* owner_space)
406      :
407      TaggedObject(),
408      owner_space_(owner_space),
409      dot_cache_(10),
410      nrm2_cache_tag_(0),
411      asum_cache_tag_(0),
412      amax_cache_tag_(0),
413      max_cache_tag_(0),
414      min_cache_tag_(0),
415      sum_cache_tag_(0),
416      sumlogs_cache_tag_(0)
417  {
418    DBG_ASSERT(IsValid(owner_space_));
419  }
420
421  inline
422  Vector* Vector::MakeNew() const
423  {
424    return owner_space_->MakeNew();
425  }
426
427  inline
428  Vector* Vector::MakeNewCopy() const
429  {
430    // ToDo: We can probably copy also the cached values for Norms etc here
431    Vector* copy = MakeNew();
432    copy->Copy(*this);
433    return copy;
434  }
435
436  inline
437  void Vector::Copy(const Vector& x)
438  {
439    CopyImpl(x);
440    ObjectChanged();
441    // Also copy any cached scalar values from the original vector
442    // ToDo: Check if that is too much overhead
443    TaggedObject::Tag x_tag = x.GetTag();
444    if (x_tag == x.nrm2_cache_tag_) {
445      nrm2_cache_tag_ = GetTag();
446      cached_nrm2_ = x.cached_nrm2_;
447    }
448    if (x_tag == x.asum_cache_tag_) {
449      asum_cache_tag_ = GetTag();
450      cached_asum_ = x.cached_asum_;
451    }
452    if (x_tag == x.amax_cache_tag_) {
453      amax_cache_tag_ = GetTag();
454      cached_amax_ = x.cached_amax_;
455    }
456    if (x_tag == x.max_cache_tag_) {
457      max_cache_tag_ = GetTag();
458      cached_max_ = x.cached_max_;
459    }
460    if (x_tag == x.min_cache_tag_) {
461      min_cache_tag_ = GetTag();
462      cached_min_ = x.cached_min_;
463    }
464    if (x_tag == x.sum_cache_tag_) {
465      sum_cache_tag_ = GetTag();
466      cached_sum_ = x.cached_sum_;
467    }
468    if (x_tag == x.sumlogs_cache_tag_) {
469      sumlogs_cache_tag_ = GetTag();
470      cached_sumlogs_ = x.cached_sumlogs_;
471    }
472  }
473
474  inline
475  void Vector::Scal(Number alpha)
476  {
477    if (alpha!=1.) {
478      TaggedObject::Tag old_tag = GetTag();
479      ScalImpl(alpha);
480      ObjectChanged();
481      if (old_tag == nrm2_cache_tag_) {
482        nrm2_cache_tag_ = GetTag();
483        cached_nrm2_ *= fabs(alpha);
484      }
485      if (old_tag == asum_cache_tag_) {
486        asum_cache_tag_ = GetTag();
487        cached_asum_ *= fabs(alpha);
488      }
489      if (old_tag == amax_cache_tag_) {
490        amax_cache_tag_ = GetTag();
491        cached_amax_ *= fabs(alpha);
492      }
493      if (old_tag == max_cache_tag_) {
494        if (alpha>=0.) {
495          max_cache_tag_ = GetTag();
496          cached_max_ *= alpha;
497        }
498        else if (alpha<0.) {
499          min_cache_tag_ = GetTag();
500          cached_min_ = cached_max_ * alpha;
501        }
502      }
503      if (old_tag == min_cache_tag_) {
504        if (alpha>=0.) {
505          min_cache_tag_ = GetTag();
506          cached_min_ *= alpha;
507        }
508        else if (alpha<0.) {
509          max_cache_tag_ = GetTag();
510          cached_max_ = cached_min_ * alpha;
511        }
512      }
513      if (old_tag == sum_cache_tag_) {
514        sum_cache_tag_ = GetTag();
515        cached_sum_ *= alpha;
516      }
517      if (old_tag == sumlogs_cache_tag_) {
518        sumlogs_cache_tag_ = GetTag();
519        cached_sumlogs_ += ((Number)Dim())*log(alpha);
520      }
521    }
522  }
523
524  inline
525  void Vector::Axpy(Number alpha, const Vector &x)
526  {
527    AxpyImpl(alpha, x);
528    ObjectChanged();
529  }
530
531  inline
532  Number Vector::Dot(const Vector &x) const
533  {
534    Number retValue;
535    if (!dot_cache_.GetCachedResult2Dep(retValue, this, &x)) {
536      retValue = DotImpl(x);
537      dot_cache_.AddCachedResult2Dep(retValue, this, &x);
538    }
539    return retValue;
540  }
541
542  inline
543  Number Vector::Nrm2() const
544  {
545    if (nrm2_cache_tag_ != GetTag()) {
546      cached_nrm2_ = Nrm2Impl();
547      nrm2_cache_tag_ = GetTag();
548    }
549    return cached_nrm2_;
550  }
551
552  inline
553  Number Vector::Asum() const
554  {
555    if (asum_cache_tag_ != GetTag()) {
556      cached_asum_ = AsumImpl();
557      asum_cache_tag_ = GetTag();
558    }
559    return cached_asum_;
560  }
561
562  inline
563  Number Vector::Amax() const
564  {
565    if (amax_cache_tag_ != GetTag()) {
566      cached_amax_ = AmaxImpl();
567      amax_cache_tag_ = GetTag();
568    }
569    return cached_amax_;
570  }
571
572  inline
573  Number Vector::Sum() const
574  {
575    if (sum_cache_tag_ != GetTag()) {
576      cached_sum_ = SumImpl();
577      sum_cache_tag_ = GetTag();
578    }
579    return cached_sum_;
580  }
581
582  inline
583  Number Vector::SumLogs() const
584  {
585    if (sumlogs_cache_tag_ != GetTag()) {
586      cached_sumlogs_ = SumLogsImpl();
587      sumlogs_cache_tag_ = GetTag();
588    }
589    return cached_sumlogs_;
590  }
591
592  inline
593  void Vector::ElementWiseSgn()
594  {
595    ElementWiseSgnImpl();
596    ObjectChanged();
597  }
598
599  inline
600  void Vector::Set(Number alpha)
601  {
602    // Could initialize caches here
603    SetImpl(alpha);
604    ObjectChanged();
605  }
606
607  inline
608  void Vector::ElementWiseDivide(const Vector& x)
609  {
610    ElementWiseDivideImpl(x);
611    ObjectChanged();
612  }
613
614  inline
615  void Vector::ElementWiseMultiply(const Vector& x)
616  {
617    ElementWiseMultiplyImpl(x);
618    ObjectChanged();
619  }
620
621  inline
622  void Vector::ElementWiseReciprocal()
623  {
624    ElementWiseReciprocalImpl();
625    ObjectChanged();
626  }
627
628  inline
629  void Vector::ElementWiseMax(const Vector& x)
630  {
631    // Could initialize some caches here
632    ElementWiseMaxImpl(x);
633    ObjectChanged();
634  }
635
636  inline
637  void Vector::ElementWiseMin(const Vector& x)
638  {
639    // Could initialize some caches here
640    ElementWiseMinImpl(x);
641    ObjectChanged();
642  }
643
644  inline
645  void Vector::ElementWiseAbs()
646  {
647    // Could initialize some caches here
648    ElementWiseAbsImpl();
649    ObjectChanged();
650  }
651
652  inline
653  void Vector::ElementWiseSqrt()
654  {
655    ElementWiseSqrtImpl();
656    ObjectChanged();
657  }
658
659  inline
660  void Vector::AddScalar(Number scalar)
661  {
662    // Could initialize some caches here
663    AddScalarImpl(scalar);
664    ObjectChanged();
665  }
666
667  inline
668  Number Vector::Max() const
669  {
670    if (max_cache_tag_ != GetTag()) {
671      cached_max_ = MaxImpl();
672      max_cache_tag_ = GetTag();
673    }
674    return cached_max_;
675  }
676
677  inline
678  Number Vector::Min() const
679  {
680    if (min_cache_tag_ != GetTag()) {
681      cached_min_ = MinImpl();
682      min_cache_tag_ = GetTag();
683    }
684    return cached_min_;
685  }
686
687  inline
688  void Vector::AddOneVector(Number a, const Vector& v1, Number c)
689  {
690    AddTwoVectorsImpl(a, v1, 0., v1, c);
691  }
692
693  inline
694  void Vector::AddTwoVectors(Number a, const Vector& v1,
695                             Number b, const Vector& v2, Number c)
696  {
697    AddTwoVectorsImpl(a, v1, b, v2, c);
698    ObjectChanged();
699  }
700
701  inline
702  Number Vector::FracToBound(const Vector& delta, Number tau) const
703  {
704    /* AW: I avoid the caching here, since it leads to overhead in the
705       quality function search.  Caches for this are in
706       CalculatedQuantities.
707    Number retValue;
708    std::vector<const TaggedObject*> tdeps(1);
709    tdeps[0] = &delta;
710    std::vector<Number> sdeps(1);
711    sdeps[0] = tau;
712    if (!frac_to_bound_cache_.GetCachedResult(retValue, tdeps, sdeps)) {
713      retValue = FracToBoundImpl(delta, tau);
714      frac_to_bound_cache_.AddCachedResult(retValue, tdeps, sdeps);
715    }
716    return retValue;
717    */
718    return FracToBoundImpl(delta, tau);
719  }
720
721  inline
722  void Vector::AddVectorQuotient(Number a, const Vector& z,
723                                 const Vector& s, Number c)
724  {
725    AddVectorQuotientImpl(a, z, s, c);
726    ObjectChanged();
727  }
728
729  inline
730  Index Vector::Dim() const
731  {
732    return owner_space_->Dim();
733  }
734
735  inline
736  SmartPtr<const VectorSpace> Vector::OwnerSpace() const
737  {
738    return owner_space_;
739  }
740
741  inline
742  VectorSpace::VectorSpace(Index dim)
743      :
744      dim_(dim)
745  {}
746
747} // namespace Ipopt
748
749// Macro definitions for debugging vectors
750#ifndef IP_DEBUG
751# define DBG_PRINT_VECTOR(__verbose_level, __vec_name, __vec)
752#else
753# define DBG_PRINT_VECTOR(__verbose_level, __vec_name, __vec) \
754   if (dbg_jrnl.Verbosity() >= (__verbose_level)) { \
755      if (dbg_jrnl.Jnlst()!=NULL) { \
756        (__vec).Print(dbg_jrnl.Jnlst(), \
757                      J_ERROR, J_DBG, \
758                      __vec_name, \
759                      dbg_jrnl.IndentationLevel()*2, \
760                      "# "); \
761      } \
762   }
763#endif // IP_DEBUG
764
765#endif
Note: See TracBrowser for help on using the repository browser.