text: // Copyright (C) 2004, 2006 International Business Machines and others.
  // All Rights Reserved.
  // This code is published under the Common Public License.
  //
  // $Id: IpMatrix.hpp 759 2006-07-07 03:07:08Z andreasw $
  //
  // Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
  
  #ifndef __IPMATRIX_HPP__
  #define __IPMATRIX_HPP__
  
  #include "IpVector.hpp"
  
  namespace SimTKIpopt
  {
  
    /* forward declarations */
    class MatrixSpace;
  
    /** Matrix Base Class. This is the base class for all derived matrix
     *  types.  All Matrices, such as Jacobian and Hessian matrices, as
     *  well as possibly the iteration matrices needed for the step
     *  computation, are of this type.
     *
     *  Deriving from Matrix:  Overload the protected XXX_Impl method.
     */
    class Matrix : public TaggedObject
    {
    public:
      /** @name Constructor/Destructor */
      //@{
      /** Constructor.  It has to be given a pointer to the
       *  corresponding MatrixSpace.
       */
      Matrix(const MatrixSpace* owner_space);
  
      /** Destructor */
      virtual ~Matrix();
      //@}
  
      /**@name Operations of the Matrix on a Vector */
      //@{
      /** Matrix-vector multiply.  Computes y = alpha * Matrix * x +
       *  beta * y.  Do not overload.  Overload MultVectorImpl instead.
       */
      void MultVector(Number alpha, const Vector& x, Number beta,
                      Vector& y) const;
  
      /** Matrix(transpose) vector multiply.  Computes y = alpha *
       *  Matrix^T * x + beta * y.  Do not overload.  Overload
       *  TransMultVectorImpl instead.
       */
      void TransMultVector(Number alpha, const Vector& x, Number beta,
                           Vector& y) const;
      //@}
  
      /** @name Methods for specialized operations.  A prototype
       *  implementation is provided, but for efficient implementation
       *  those should be specially implemented.
       */
      //@{
      /** X = X + alpha*(Matrix S^{-1} Z).  Should be implemented
       *  efficiently for the ExansionMatrix
       */
      void AddMSinvZ(Number alpha, const Vector& S, const Vector& Z,
                     Vector& X) const;
  
      /** X = S^{-1} (r + alpha*Z*M^Td).   Should be implemented
       *  efficiently for the ExansionMatrix
       */
      void SinvBlrmZMTdBr(Number alpha, const Vector& S,
                          const Vector& R, const Vector& Z,
                          const Vector& D, Vector& X) const;
      //@}
  
      /** Method for determining if all stored numbers are valid (i.e.,
       *  no Inf or Nan). */
      bool HasValidNumbers() const;
  
      //* @name Information about the size of the matrix */
      //@{
      /** Number of rows */
      Index  NRows() const;
  
      /** Number of columns */
      Index  NCols() const;
      //@}
  
      /** Print detailed information about the matrix. Do not overload.
       *  Overload PrintImpl instead.
       */
      //@{
      virtual void Print(SmartPtr<const Journalist> jnlst,
                         EJournalLevel level,
                         EJournalCategory category,
                         const std::string& name,
                         Index indent=0,
                         const std::string& prefix="") const;
      virtual void Print(const Journalist& jnlst,
                         EJournalLevel level,
                         EJournalCategory category,
                         const std::string& name,
                         Index indent=0,
                         const std::string& prefix="") const;
      //@}
  
      /** Return the owner MatrixSpace*/
      SmartPtr<const MatrixSpace> OwnerSpace() const;
  
    protected:
      /** @name implementation methods (derived classes MUST
       *  overload these pure virtual protected methods.
       */
      //@{
      /** Matrix-vector multiply.  Computes y = alpha * Matrix * x +
       *  beta * y
       */
      virtual void MultVectorImpl(Number alpha, const Vector& x, Number beta, Vector& y) const =0;
  
      /** Matrix(transpose) vector multiply.
       * Computes y = alpha * Matrix^T * x  +  beta * y
       */
      virtual void TransMultVectorImpl(Number alpha, const Vector& x, Number beta, Vector& y) const =0;
  
      /** X = X + alpha*(Matrix S^{-1} Z).  Prototype for this
       *  specialize method is provided, but for efficient
       *  implementation it should be overloaded for the expansion matrix.
       */
      virtual void AddMSinvZImpl(Number alpha, const Vector& S, const Vector& Z,
                                 Vector& X) const;
  
      /** X = S^{-1} (r + alpha*Z*M^Td).   Should be implemented
       *  efficiently for the ExpansionMatrix.
       */
      virtual void SinvBlrmZMTdBrImpl(Number alpha, const Vector& S,
                                      const Vector& R, const Vector& Z,
                                      const Vector& D, Vector& X) const;
  
      /** Method for determining if all stored numbers are valid (i.e.,
       *  no Inf or Nan). A default implementation always returning true
       *  is provided, but if possible it should be implemented. */
      virtual bool HasValidNumbersImpl() const
      {
        return true;
      }
  
      /** Print detailed information about the matrix. */
      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 */
      Matrix();
  
      /** Copy constructor */
      Matrix(const Matrix&);
  
      /** Overloaded Equals Operator */
      Matrix& operator=(const Matrix&);
      //@}
  
      const SmartPtr<const MatrixSpace> owner_space_;
  
      /**@name CachedResults data members */
      //@{
      mutable TaggedObject::Tag valid_cache_tag_;
      mutable bool cached_valid_;
      //@}
    };
  
  
    /** MatrixSpace base class, corresponding to the Matrix base class.
     *  For each Matrix implementation, a corresponding MatrixSpace has
     *  to be implemented.  A MatrixSpace is able to create new Matrices
     *  of a specific type.  The MatrixSpace should also store
     *  information that is common to all Matrices of that type.  For
     *  example, the dimensions of a Matrix is stored in the MatrixSpace
     *  base class.
     */
    class MatrixSpace : public ReferencedObject
    {
    public:
      /** @name Constructors/Destructors */
      //@{
      /** Constructor, given the number rows and columns of all matrices
       *  generated by this MatrixSpace.
       */
      MatrixSpace(Index nRows, Index nCols);
  
      /** Destructor */
      virtual ~MatrixSpace();
      //@}
  
      /** Pure virtual method for creating a new Matrix of the
       *  corresponding type.
       */
      virtual Matrix* MakeNew() const=0;
  
      /** Accessor function for the number of rows. */
      Index NRows() const
      {
        return nRows_;
      }
      /** Accessor function for the number of columns. */
      Index NCols() const
      {
        return nCols_;
      }
  
      /** Method to test if a given matrix belongs to a particular
       *  matrix space.
       */
      bool IsMatrixFromSpace(const Matrix& matrix) const;
  
    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 */
      MatrixSpace();
  
      /** Copy constructor */
      MatrixSpace(const MatrixSpace&);
  
      /** Overloaded Equals Operator */
      MatrixSpace& operator=(const MatrixSpace&);
      //@}
  
      /** Number of rows for all matrices of this type. */
      const Index nRows_;
      /** Number of columns for all matrices of this type. */
      const Index nCols_;
    };
  
  
    /* Inline Methods */
    inline
    Matrix::Matrix(const MatrixSpace* owner_space)
        :
        TaggedObject(),
        owner_space_(owner_space),
        valid_cache_tag_(0)
    {}
  
    inline
    Matrix::~Matrix()
    {}
  
    inline
    Index  Matrix::NRows() const
    {
      return owner_space_->NRows();
    }
  
    inline
    Index  Matrix::NCols() const
    {
      return owner_space_->NCols();
    }
  
    inline
    SmartPtr<const MatrixSpace> Matrix::OwnerSpace() const
    {
      return owner_space_;
    }
  
  } // namespace Ipopt
  
  // Macro definitions for debugging matrices
  #ifndef IP_DEBUG
  # define DBG_PRINT_MATRIX(__verbose_level, __mat_name, __mat)
  #else
  # define DBG_PRINT_MATRIX(__verbose_level, __mat_name, __mat) \
     if (dbg_jrnl.Verbosity() >= (__verbose_level)) { \
        if (dbg_jrnl.Jnlst()!=NULL) { \
          (__mat).Print(dbg_jrnl.Jnlst(), \
                J_ERROR, J_DBG, \
                __mat_name, \
                dbg_jrnl.IndentationLevel()*2, \
                "# "); \
        } \
     }
  #endif // IP_DEBUG
  
  #endif
  