// 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 __IPMATRIX_HPP__
#define __IPMATRIX_HPP__

#include "IpVector.hpp"

namespace Ipopt
{

/* 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 IPOPTLIB_EXPORT Matrix: public TaggedObject
{
public:
   /** @name Constructor/Destructor */
   ///@{
   /** Constructor.
    *
    * It has to be given a pointer to the corresponding MatrixSpace.
    */
   Matrix(
      const MatrixSpace* owner_space
   )
      : TaggedObject(),
        owner_space_(owner_space),
        valid_cache_tag_(0)
   { }

   /** Destructor */
   virtual ~Matrix()
   { }
   ///@}

   /**@name Operations of the Matrix on a Vector */
   ///@{
   /** Matrix-vector multiply.
    *
    * Computes y = alpha * Matrix * x + beta * y.
    *
    * @attention Do not overload. Overload MultVectorImpl instead.
    */
   void MultVector(
      Number        alpha,
      const Vector& x,
      Number        beta,
      Vector&       y
   ) const
   {
      MultVectorImpl(alpha, x, beta, y);
   }

   /** Matrix(transpose) vector multiply.
    *
    *  Computes y = alpha * Matrix^T * x + beta * y.
    *
    *  @attention Do not overload.  Overload TransMultVectorImpl instead.
    */
   void TransMultVector(
      Number        alpha,
      const Vector& x,
      Number        beta,
      Vector&       y
   ) const
   {
      TransMultVectorImpl(alpha, x, beta, y);
   }
   ///@}

   /** @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 ExpansionMatrix.
    */
   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 ExpansionMatrix.
    */
   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 */
   inline Index NRows() const;

   /** Number of columns */
   inline Index NCols() const;
   ///@}

   /** @name Norms of the individual rows and columns */
   ///@{
   /** Compute the max-norm of the rows in the matrix.
    *
    *  The result is stored in rows_norms.
    *  The vector is assumed to be initialized if init is false.
    */
   void ComputeRowAMax(
      Vector& rows_norms,
      bool    init = true
   ) const
   {
      DBG_ASSERT(NRows() == rows_norms.Dim());
      if( init )
      {
         rows_norms.Set(0.);
      }
      ComputeRowAMaxImpl(rows_norms, init);
   }

   /** Compute the max-norm of the columns in the matrix.
    *
    *  The result is stored in cols_norms.
    *  The vector is assumed to be initialized if init is false.
    */
   void ComputeColAMax(
      Vector& cols_norms,
      bool    init = true
   ) const
   {
      DBG_ASSERT(NCols() == cols_norms.Dim());
      if( init )
      {
         cols_norms.Set(0.);
      }
      ComputeColAMaxImpl(cols_norms, init);
   }
   ///@}

   /** Print detailed information about the matrix.
    *
    * @attention 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 */
   inline 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;
   }

   /** Compute the max-norm of the rows in the matrix.
    *
    *  The result is stored in rows_norms.
    *  The vector is assumed to be initialized if init is false. */
   virtual void ComputeRowAMaxImpl(
      Vector& rows_norms,
      bool    init
   ) const = 0;

   /** Compute the max-norm of the columns in the matrix.
    *
    *  The result is stored in cols_norms.
    *  The vector is assumed to be initialized if init is false. */
   virtual void ComputeColAMaxImpl(
      Vector& cols_norms,
      bool    init
   ) const = 0;

   /** 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&
   );

   /** Default Assignment 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 IPOPTLIB_EXPORT 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)
      : nRows_(nRows),
        nCols_(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
   {
      return (matrix.OwnerSpace() == this);
   }

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&
   );

   /** Default Assignment 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 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
#if IPOPT_VERBOSITY == 0
# 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 // #if IPOPT_VERBOSITY == 0

#endif