// Copyright (C) 2004, 2008 International Business Machines and others. // All Rights Reserved. // This code is published under the Common Public License. // // $Id: IpVector.hpp 1327 2008-09-18 19:01:17Z andreasw $ // // Authors: Carl Laird, Andreas Waechter IBM 2004-08-13 #ifndef __IPVECTOR_HPP__ #define __IPVECTOR_HPP__ #include "IpTypes.hpp" #include "IpTaggedObject.hpp" #include "IpCachedResults.hpp" #include "IpSmartPtr.hpp" #include "IpJournalist.hpp" #include "IpException.hpp" #include namespace Ipopt { /** Exception that can be used to flag unimplemented linear algebra * methods */ DECLARE_STD_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED); /* forward declarations */ class VectorSpace; /** Vector Base Class. * This is the base class for all derived vector types. Those vectors * are meant to store entities like iterates, Lagrangian multipliers, * constraint values etc. The implementation of a vector type depends * on the computational environment (e.g. just a double array on a shared * memory machine, or distributed double arrays for a distributed * memory machine.) * * Deriving from Vector: This class inherits from tagged object to * implement an advanced caching scheme. Because of this, the * TaggedObject method ObjectChanged() must be called each time the * Vector changes. If you overload the XXXX_Impl protected methods, * this taken care of (along with caching if possible) for you. If * you have additional methods in your derived class that change the * underlying data (vector values), you MUST remember to call * ObjectChanged() AFTER making the change! */ class Vector : public TaggedObject { public: /** @name Constructor/Destructor */ //@{ /** Constructor. It has to be given a pointer to the * corresponding VectorSpace. */ Vector(const VectorSpace* owner_space); /** Destructor */ virtual ~Vector(); //@} /** Create new Vector of the same type with uninitialized data */ Vector* MakeNew() const; /** Create new Vector of the same type and copy the data over */ Vector* MakeNewCopy() const; /**@name Standard BLAS-1 Operations * (derived classes do NOT overload these * methods, instead, overload the * protected versions of these methods). */ //@{ /** Copy the data of the vector x into this vector (DCOPY). */ void Copy(const Vector& x); /** Scales the vector by scalar alpha (DSCAL) */ void Scal(Number alpha); /** Add the multiple alpha of vector x to this vector (DAXPY) */ void Axpy(Number alpha, const Vector &x); /** Computes inner product of vector x with this (DDOT) */ Number Dot(const Vector &x) const; /** Computes the 2-norm of this vector (DNRM2) */ Number Nrm2() const; /** Computes the 1-norm of this vector (DASUM) */ Number Asum() const; /** Computes the max-norm of this vector (based on IDAMAX) */ Number Amax() const; //@} /** @name Additional (Non-BLAS) Vector Methods * (derived classes do NOT overload these * methods, instead, overload the * protected versions of these methods). */ //@{ /** Set each element in the vector to the scalar alpha. */ void Set(Number alpha); /** Element-wise division \f$y_i \gets y_i/x_i\f$*/ void ElementWiseDivide(const Vector& x); /** Element-wise multiplication \f$y_i \gets y_i*x_i\f$ */ void ElementWiseMultiply(const Vector& x); /** Element-wise max against entries in x */ void ElementWiseMax(const Vector& x); /** Element-wise min against entries in x */ void ElementWiseMin(const Vector& x); /** Reciprocates the entries in the vector */ void ElementWiseReciprocal(); /** Absolute values of the entries in the vector */ void ElementWiseAbs(); /** Element-wise square root of the entries in the vector */ void ElementWiseSqrt(); /** Replaces the vector values with their sgn values ( -1 if x_i < 0, 0 if x_i == 0, and 1 if x_i > 0) */ void ElementWiseSgn(); /** Add scalar to every vector component */ void AddScalar(Number scalar); /** Returns the maximum value in the vector */ Number Max() const; /** Returns the minimum value in the vector */ Number Min() const; /** Returns the sum of the vector entries */ Number Sum() const; /** Returns the sum of the logs of each vector entry */ Number SumLogs() const; //@} /** @name Methods for specialized operations. A prototype * implementation is provided, but for efficient implementation * those should be specially implemented. */ //@{ /** Add one vector, y = a * v1 + c * y. This is automatically * reduced to call AddTwoVectors. */ void AddOneVector(Number a, const Vector& v1, Number c); /** Add two vectors, y = a * v1 + b * v2 + c * y. Here, this * vector is y */ void AddTwoVectors(Number a, const Vector& v1, Number b, const Vector& v2, Number c); /** Fraction to the boundary parameter. Computes \f$\alpha = * \max\{\bar\alpha\in(0,1] : x + \bar\alpha \Delta \geq (1-\tau)x\}\f$ */ Number FracToBound(const Vector& delta, Number tau) const; /** Add the quotient of two vectors, y = a * z/s + c * y. */ void AddVectorQuotient(Number a, const Vector& z, const Vector& s, Number c); //@} /** Method for determining if all stored numbers are valid (i.e., * no Inf or Nan). */ bool HasValidNumbers() const; /** @name Accessor methods */ //@{ /** Dimension of the Vector */ Index Dim() const; /** Return the owner VectorSpace*/ SmartPtr OwnerSpace() const; //@} /** @name Output methods * (derived classes do NOT overload these * methods, instead, overload the * protected versions of these methods). */ //@{ /** Print the entire vector */ void Print(SmartPtr jnlst, EJournalLevel level, EJournalCategory category, const std::string& name, Index indent=0, const std::string& prefix="") const; void Print(const Journalist& jnlst, EJournalLevel level, EJournalCategory category, const std::string& name, Index indent=0, const std::string& prefix="") const; //@} protected: /** @name implementation methods (derived classes MUST * overload these pure virtual protected methods.) */ //@{ /** Copy the data of the vector x into this vector (DCOPY). */ virtual void CopyImpl(const Vector& x)=0; /** Scales the vector by scalar alpha (DSCAL) */ virtual void ScalImpl(Number alpha)=0; /** Add the multiple alpha of vector x to this vector (DAXPY) */ virtual void AxpyImpl(Number alpha, const Vector &x)=0; /** Computes inner product of vector x with this (DDOT) */ virtual Number DotImpl(const Vector &x) const =0; /** Computes the 2-norm of this vector (DNRM2) */ virtual Number Nrm2Impl() const =0; /** Computes the 1-norm of this vector (DASUM) */ virtual Number AsumImpl() const =0; /** Computes the max-norm of this vector (based on IDAMAX) */ virtual Number AmaxImpl() const =0; /** Set each element in the vector to the scalar alpha. */ virtual void SetImpl(Number alpha)=0; /** Element-wise division \f$y_i \gets y_i/x_i\f$*/ virtual void ElementWiseDivideImpl(const Vector& x)=0; /** Element-wise multiplication \f$y_i \gets y_i*x_i\f$ */ virtual void ElementWiseMultiplyImpl(const Vector& x)=0; /** Element-wise max against entries in x */ virtual void ElementWiseMaxImpl(const Vector& x)=0; /** Element-wise min against entries in x */ virtual void ElementWiseMinImpl(const Vector& x)=0; /** Reciprocates the elements of the vector */ virtual void ElementWiseReciprocalImpl()=0; /** Take elementwise absolute values of the elements of the vector */ virtual void ElementWiseAbsImpl()=0; /** Take elementwise square-root of the elements of the vector */ virtual void ElementWiseSqrtImpl()=0; /** Replaces entries with sgn of the entry */ virtual void ElementWiseSgnImpl()=0; /** Add scalar to every component of vector */ virtual void AddScalarImpl(Number scalar)=0; /** Max value in the vector */ virtual Number MaxImpl() const=0; /** Min number in the vector */ virtual Number MinImpl() const=0; /** Sum of entries in the vector */ virtual Number SumImpl() const=0; /** Sum of logs of entries in the vector */ virtual Number SumLogsImpl() const=0; /** Add two vectors (a * v1 + b * v2). Result is stored in this vector. */ virtual void AddTwoVectorsImpl(Number a, const Vector& v1, Number b, const Vector& v2, Number c); /** Fraction to boundary parameter. */ virtual Number FracToBoundImpl(const Vector& delta, Number tau) const; /** Add the quotient of two vectors */ virtual void AddVectorQuotientImpl(Number a, const Vector& z, const Vector& s, Number c); /** Method for determining if all stored numbers are valid (i.e., * no Inf or Nan). A default implementation using Asum is * provided. */ virtual bool HasValidNumbersImpl() const; /** Print the entire vector */ virtual void PrintImpl(const Journalist& jnlst, EJournalLevel level, EJournalCategory category, const std::string& name, Index indent, const std::string& prefix) const =0; //@} private: /**@name Default Compiler Generated Methods * (Hidden to avoid implicit creation/calling). * These methods are not implemented and * we do not want the compiler to implement * them for us, so we declare them private * and do not define them. This ensures that * they will not be implicitly created/called. */ //@{ /** Default constructor */ Vector(); /** Copy constructor */ Vector(const Vector&); /** Overloaded Equals Operator */ Vector& operator=(const Vector&); //@} /** Vector Space */ const SmartPtr owner_space_; /**@name CachedResults data members */ //@{ /** Cache for dot products */ mutable CachedResults dot_cache_; mutable TaggedObject::Tag nrm2_cache_tag_; mutable Number cached_nrm2_; mutable TaggedObject::Tag asum_cache_tag_; mutable Number cached_asum_; mutable TaggedObject::Tag amax_cache_tag_; mutable Number cached_amax_; mutable TaggedObject::Tag max_cache_tag_; mutable Number cached_max_; mutable TaggedObject::Tag min_cache_tag_; mutable Number cached_min_; mutable TaggedObject::Tag sum_cache_tag_; mutable Number cached_sum_; mutable TaggedObject::Tag sumlogs_cache_tag_; mutable Number cached_sumlogs_; mutable TaggedObject::Tag valid_cache_tag_; mutable bool cached_valid_; // AW: I removed this cache since it gets in the way for the // quality function search // /** Cache for FracToBound */ // mutable CachedResults frac_to_bound_cache_; //@} }; /** VectorSpace base class, corresponding to the Vector base class. * For each Vector implementation, a corresponding VectorSpace has * to be implemented. A VectorSpace is able to create new Vectors * of a specific type. The VectorSpace should also store * information that is common to all Vectors of that type. For * example, the dimension of a Vector is stored in the VectorSpace * base class. */ class VectorSpace : public ReferencedObject { public: /** @name Constructors/Destructors */ //@{ /** Constructor, given the dimension of all vectors generated by * this VectorSpace. */ VectorSpace(Index dim); /** Destructor */ virtual ~VectorSpace() {} //@} /** Pure virtual method for creating a new Vector of the * corresponding type. */ virtual Vector* MakeNew() const=0; /** Accessor function for the dimension of the vectors of this type.*/ Index Dim() const { return dim_; } private: /**@name Default Compiler Generated Methods * (Hidden to avoid implicit creation/calling). * These methods are not implemented and * we do not want the compiler to implement * them for us, so we declare them private * and do not define them. This ensures that * they will not be implicitly created/called. */ //@{ /** default constructor */ VectorSpace(); /** Copy constructor */ VectorSpace(const VectorSpace&); /** Overloaded Equals Operator */ VectorSpace& operator=(const VectorSpace&); //@} /** Dimension of the vectors in this vector space. */ const Index dim_; }; /* inline methods */ inline Vector::~Vector() {} inline Vector::Vector(const VectorSpace* owner_space) : TaggedObject(), owner_space_(owner_space), dot_cache_(10), nrm2_cache_tag_(0), asum_cache_tag_(0), amax_cache_tag_(0), max_cache_tag_(0), min_cache_tag_(0), sum_cache_tag_(0), sumlogs_cache_tag_(0), cached_valid_(0) { DBG_ASSERT(IsValid(owner_space_)); } inline Vector* Vector::MakeNew() const { return owner_space_->MakeNew(); } inline Vector* Vector::MakeNewCopy() const { // ToDo: We can probably copy also the cached values for Norms etc here Vector* copy = MakeNew(); copy->Copy(*this); return copy; } inline void Vector::Copy(const Vector& x) { CopyImpl(x); ObjectChanged(); // Also copy any cached scalar values from the original vector // ToDo: Check if that is too much overhead TaggedObject::Tag x_tag = x.GetTag(); if (x_tag == x.nrm2_cache_tag_) { nrm2_cache_tag_ = GetTag(); cached_nrm2_ = x.cached_nrm2_; } if (x_tag == x.asum_cache_tag_) { asum_cache_tag_ = GetTag(); cached_asum_ = x.cached_asum_; } if (x_tag == x.amax_cache_tag_) { amax_cache_tag_ = GetTag(); cached_amax_ = x.cached_amax_; } if (x_tag == x.max_cache_tag_) { max_cache_tag_ = GetTag(); cached_max_ = x.cached_max_; } if (x_tag == x.min_cache_tag_) { min_cache_tag_ = GetTag(); cached_min_ = x.cached_min_; } if (x_tag == x.sum_cache_tag_) { sum_cache_tag_ = GetTag(); cached_sum_ = x.cached_sum_; } if (x_tag == x.sumlogs_cache_tag_) { sumlogs_cache_tag_ = GetTag(); cached_sumlogs_ = x.cached_sumlogs_; } } inline void Vector::Axpy(Number alpha, const Vector &x) { AxpyImpl(alpha, x); ObjectChanged(); } inline Number Vector::Dot(const Vector &x) const { // The current implementation of the caching doesn't allow to have // a dependency of something with itself. Therefore, we use the // Nrm2 method if the dot product is to be taken with the vector // itself. Might be more efficient anyway. if (this==&x) { Number nrm2 = Nrm2(); return nrm2*nrm2; } Number retValue; if (!dot_cache_.GetCachedResult2Dep(retValue, this, &x)) { retValue = DotImpl(x); dot_cache_.AddCachedResult2Dep(retValue, this, &x); } return retValue; } inline Number Vector::Nrm2() const { if (nrm2_cache_tag_ != GetTag()) { cached_nrm2_ = Nrm2Impl(); nrm2_cache_tag_ = GetTag(); } return cached_nrm2_; } inline Number Vector::Asum() const { if (asum_cache_tag_ != GetTag()) { cached_asum_ = AsumImpl(); asum_cache_tag_ = GetTag(); } return cached_asum_; } inline Number Vector::Amax() const { if (amax_cache_tag_ != GetTag()) { cached_amax_ = AmaxImpl(); amax_cache_tag_ = GetTag(); } return cached_amax_; } inline Number Vector::Sum() const { if (sum_cache_tag_ != GetTag()) { cached_sum_ = SumImpl(); sum_cache_tag_ = GetTag(); } return cached_sum_; } inline Number Vector::SumLogs() const { if (sumlogs_cache_tag_ != GetTag()) { cached_sumlogs_ = SumLogsImpl(); sumlogs_cache_tag_ = GetTag(); } return cached_sumlogs_; } inline void Vector::ElementWiseSgn() { ElementWiseSgnImpl(); ObjectChanged(); } inline void Vector::Set(Number alpha) { // Could initialize caches here SetImpl(alpha); ObjectChanged(); } inline void Vector::ElementWiseDivide(const Vector& x) { ElementWiseDivideImpl(x); ObjectChanged(); } inline void Vector::ElementWiseMultiply(const Vector& x) { ElementWiseMultiplyImpl(x); ObjectChanged(); } inline void Vector::ElementWiseReciprocal() { ElementWiseReciprocalImpl(); ObjectChanged(); } inline void Vector::ElementWiseMax(const Vector& x) { // Could initialize some caches here ElementWiseMaxImpl(x); ObjectChanged(); } inline void Vector::ElementWiseMin(const Vector& x) { // Could initialize some caches here ElementWiseMinImpl(x); ObjectChanged(); } inline void Vector::ElementWiseAbs() { // Could initialize some caches here ElementWiseAbsImpl(); ObjectChanged(); } inline void Vector::ElementWiseSqrt() { ElementWiseSqrtImpl(); ObjectChanged(); } inline void Vector::AddScalar(Number scalar) { // Could initialize some caches here AddScalarImpl(scalar); ObjectChanged(); } inline Number Vector::Max() const { if (max_cache_tag_ != GetTag()) { cached_max_ = MaxImpl(); max_cache_tag_ = GetTag(); } return cached_max_; } inline Number Vector::Min() const { if (min_cache_tag_ != GetTag()) { cached_min_ = MinImpl(); min_cache_tag_ = GetTag(); } return cached_min_; } inline void Vector::AddOneVector(Number a, const Vector& v1, Number c) { AddTwoVectors(a, v1, 0., v1, c); } inline void Vector::AddTwoVectors(Number a, const Vector& v1, Number b, const Vector& v2, Number c) { AddTwoVectorsImpl(a, v1, b, v2, c); ObjectChanged(); } inline Number Vector::FracToBound(const Vector& delta, Number tau) const { /* AW: I avoid the caching here, since it leads to overhead in the quality function search. Caches for this are in CalculatedQuantities. Number retValue; std::vector tdeps(1); tdeps[0] = δ std::vector sdeps(1); sdeps[0] = tau; if (!frac_to_bound_cache_.GetCachedResult(retValue, tdeps, sdeps)) { retValue = FracToBoundImpl(delta, tau); frac_to_bound_cache_.AddCachedResult(retValue, tdeps, sdeps); } return retValue; */ return FracToBoundImpl(delta, tau); } inline void Vector::AddVectorQuotient(Number a, const Vector& z, const Vector& s, Number c) { AddVectorQuotientImpl(a, z, s, c); ObjectChanged(); } inline bool Vector::HasValidNumbers() const { if (valid_cache_tag_ != GetTag()) { cached_valid_ = HasValidNumbersImpl(); valid_cache_tag_ = GetTag(); } return cached_valid_; } inline Index Vector::Dim() const { return owner_space_->Dim(); } inline SmartPtr Vector::OwnerSpace() const { return owner_space_; } inline VectorSpace::VectorSpace(Index dim) : dim_(dim) {} } // namespace Ipopt // Macro definitions for debugging vectors #if COIN_IPOPT_VERBOSITY == 0 # define DBG_PRINT_VECTOR(__verbose_level, __vec_name, __vec) #else # define DBG_PRINT_VECTOR(__verbose_level, __vec_name, __vec) \ if (dbg_jrnl.Verbosity() >= (__verbose_level)) { \ if (dbg_jrnl.Jnlst()!=NULL) { \ (__vec).Print(dbg_jrnl.Jnlst(), \ J_ERROR, J_DBG, \ __vec_name, \ dbg_jrnl.IndentationLevel()*2, \ "# "); \ } \ } #endif //if COIN_IPOPT_VERBOSITY == 0 #endif