ipopt-src 0.2.3+3.14.16

// Copyright (C) 2004, 2008 International Business Machines and others.
// All Rights Reserved.
// This code is published under the Eclipse Public License.
//
// 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 <vector>

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 IPOPTLIB_EXPORT Vector: public TaggedObject
{
public:
   /** @name Constructor/Destructor */
   ///@{
   /** Constructor.
    *
    *  It has to be given a pointer to the corresponding VectorSpace.
    */
   inline Vector(
      const VectorSpace* owner_space
   );

   /** Destructor */
   inline virtual ~Vector();
   ///@}

   /** Create new Vector of the same type with uninitialized data */
   inline Vector* MakeNew() const;

   /** Create new Vector of the same type and copy the data over */
   inline 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). */
   inline 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) */
   inline void Axpy(
      Number        alpha,
      const Vector& x
   );

   /** Computes inner product of vector x with this (DDOT) */
   inline Number Dot(
      const Vector& x
   ) const;

   /** Computes the 2-norm of this vector (DNRM2) */
   inline Number Nrm2() const;

   /** Computes the 1-norm of this vector (DASUM) */
   inline Number Asum() const;

   /** Computes the max-norm of this vector (based on IDAMAX) */
   inline 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. */
   inline void Set(
      Number alpha
   );

   /** Element-wise division  \f$y_i \gets y_i/x_i\f$ */
   inline void ElementWiseDivide(
      const Vector& x
   );

   /** Element-wise multiplication \f$y_i \gets y_i*x_i\f$ */
   inline void ElementWiseMultiply(
      const Vector& x
   );

   /** Element-wise selection \f$y_i \gets sgn(y_i)x_i, if y_i != 0, else 0\f$
    *
    * So like ElementWiseSgn().ElementWiseMultiply(x), but with 0*inf = 0.
    * @since 3.14.0
    */
   inline void ElementWiseSelect(
      const Vector& x
   );

   /** Element-wise max against entries in x */
   inline void ElementWiseMax(
      const Vector& x
   );

   /** Element-wise min against entries in x */
   inline void ElementWiseMin(
      const Vector& x
   );

   /** Reciprocates the entries in the vector */
   inline void ElementWiseReciprocal();

   /** Absolute values of the entries in the vector */
   inline void ElementWiseAbs();

   /** Element-wise square root of the entries in the vector */
   inline 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)
    */
   inline void ElementWiseSgn();

   /** Add scalar to every vector component */
   inline void AddScalar(
      Number scalar
   );

   /** Returns the maximum value in the vector */
   inline Number Max() const;

   /** Returns the minimum value in the vector */
   inline Number Min() const;

   /** Returns the sum of the vector entries */
   inline Number Sum() const;

   /** Returns the sum of the logs of each vector entry */
   inline 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.
    */
   inline void AddOneVector(
      Number        a,
      const Vector& v1,
      Number        c
   );

   /** Add two vectors, y = a * v1 + b * v2 + c * y.
    *
    * Here, this vector is y.
    */
   inline 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$
    */
   inline Number FracToBound(
      const Vector& delta,
      Number        tau
   ) const;

   /** Add the quotient of two vectors, y = a * z/s + c * y. */
   inline 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). */
   inline bool HasValidNumbers() const;

   /** @name Accessor methods */
   ///@{
   /** Dimension of the Vector */
   inline Index Dim() const;

   /** Return the owner VectorSpace*/
   inline SmartPtr<const VectorSpace> 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<const Journalist> 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 selection \f$y_i \gets sgn(y_i)x_i, if y_i != 0, else 0\f$ */
   virtual void ElementWiseSelectImpl(
      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&
   );

   /** Default Assignment Operator */
   Vector& operator=(
      const Vector&
   );
   ///@}

   /** Vector Space */
   const SmartPtr<const VectorSpace> owner_space_;

   /**@name CachedResults data members */
   ///@{
   /** Cache for dot products */
   mutable CachedResults<Number> 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<Number> 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 IPOPTLIB_EXPORT 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&
   );

   /** Default Assignment 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::ElementWiseSelect(
   const Vector& x
)
{
   ElementWiseSelectImpl(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<const TaggedObject*> tdeps(1);
    tdeps[0] = &delta;
    std::vector<Number> 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<const VectorSpace> Vector::OwnerSpace() const
{
   return owner_space_;
}

inline VectorSpace::VectorSpace(
   Index dim
)
   : dim_(dim)
{ }

} // namespace Ipopt

// Macro definitions for debugging vectors
#if 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 IPOPT_VERBOSITY == 0

#endif