source: branches/dev/LinAlg/IpVector.hpp @ 501

Last change on this file since 501 was 501, checked in by claird, 15 years ago

Cleaned up all the Copyright comments.

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