FGMRES refinement for Primal Dual System

src/Algorithm/IpPDFullSpaceSolver.hpp

@@ -132,6 +132,9 @@ class PDFullSpaceSolver: public PDSystemSolver
     */
    Number residual_improvement_factor_;
 
+   /** Use (flexible) GMRES refinement instead of iterative refinement */
+   bool gmres_refinement_;
+
    /** Tolerance for heuristic to ignore wrong inertia */
    Number neg_curv_test_tol_;
 
@@ -208,7 +211,8 @@ class PDFullSpaceSolver: public PDSystemSolver
    Number ComputeResidualRatio(
       const IteratesVector& rhs,
       const IteratesVector& res,
-      const IteratesVector& resid
+      const IteratesVector& resid,
+      Number AInfNrm
    );
 
    /** @name Auxiliary functions */
@@ -224,6 +228,62 @@ class PDFullSpaceSolver: public PDSystemSolver
       Vector&       X
    );
    ///@}
+
+   bool SolveGMRES(
+      Number                alpha,
+      Number                beta,
+      const IteratesVector& rhs,
+      IteratesVector&       res,
+      bool                  allow_inexact = false,
+      bool                  improve_solution = false
+   );
+
+   bool GMRES(
+      const SymMatrix&      W,
+      const Matrix&         J_c,
+      const Matrix&         J_d,
+      const Matrix&         Px_L,
+      const Matrix&         Px_U,
+      const Matrix&         Pd_L,
+      const Matrix&         Pd_U,
+      const Vector&         z_L,
+      const Vector&         z_U,
+      const Vector&         v_L,
+      const Vector&         v_U,
+      const Vector&         slack_x_L,
+      const Vector&         slack_x_U,
+      const Vector&         slack_s_L,
+      const Vector&         slack_s_U,
+      const Vector&         sigma_x,
+      const Vector&         sigma_s,
+      const IteratesVector& rhs,
+      IteratesVector&       res,
+      IteratesVector&       resid,
+      bool                  improve_solution,
+      bool                  resolve_with_better_quality,
+      bool                  pretend_singular,
+      bool&                 solve_retval
+   );
+
+   Number NrmInf(
+      const SymMatrix& W,
+      const Matrix&    J_c,
+      const Matrix&    J_d,
+      const Matrix&    Px_L,
+      const Matrix&    Px_U,
+      const Matrix&    Pd_L,
+      const Matrix&    Pd_U,
+      const Vector&    z_L,
+      const Vector&    z_U,
+      const Vector&    v_L,
+      const Vector&    v_U,
+      const Vector&    slack_x_L,
+      const Vector&    slack_x_U,
+      const Vector&    slack_s_L,
+      const Vector&    slack_s_U,
+      IteratesVector&  tmp
+   );
+
 };
 
 } // namespace Ipopt

src/LinAlg/IpCompoundMatrix.cpp

@@ -715,6 +715,118 @@ void CompoundMatrix::ComputeColAMaxImpl(
    }
 }
 
+void CompoundMatrix::ComputeRowA1Impl(
+   Vector& rows_norms,
+   bool    /*init*/
+) const
+{
+   if( !matrices_valid_ )
+   {
+      matrices_valid_ = MatricesValid();
+   }
+   DBG_ASSERT(matrices_valid_);
+
+   // The vector is assumed to be compound Vectors as well except if
+   // there is only one component
+   CompoundVector* comp_vec = dynamic_cast<CompoundVector*>(&rows_norms);
+
+#ifndef ALLOW_NESTED
+   //  A few sanity checks
+   if (comp_vec)
+   {
+      DBG_ASSERT(NComps_Rows() == comp_vec->NComps());
+   }
+   else
+   {
+      DBG_ASSERT(NComps_Rows() == 1);
+   }
+#endif
+   if( comp_vec )
+   {
+      if( NComps_Rows() != comp_vec->NComps() )
+      {
+         comp_vec = NULL;
+      }
+   }
+
+   for( Index jcol = 0; jcol < NComps_Cols(); jcol++ )
+   {
+      for( Index irow = 0; irow < NComps_Rows(); irow++ )
+      {
+         if( ConstComp(irow, jcol) )
+         {
+            SmartPtr<Vector> vec_i;
+            if( comp_vec )
+            {
+               vec_i = comp_vec->GetCompNonConst(irow);
+            }
+            else
+            {
+               vec_i = &rows_norms;
+            }
+            DBG_ASSERT(IsValid(vec_i));
+            ConstComp(irow, jcol)->ComputeRowA1(*vec_i, false);
+         }
+      }
+   }
+}
+
+void CompoundMatrix::ComputeColA1Impl(
+   Vector& cols_norms,
+   bool    /*init*/
+) const
+{
+   if( !matrices_valid_ )
+   {
+      matrices_valid_ = MatricesValid();
+   }
+   DBG_ASSERT(matrices_valid_);
+
+   // The vector is assumed to be compound Vectors as well except if
+   // there is only one component
+   CompoundVector* comp_vec = dynamic_cast<CompoundVector*>(&cols_norms);
+
+#ifndef ALLOW_NESTED
+   //  A few sanity checks
+   if (comp_vec)
+   {
+      DBG_ASSERT(NComps_Cols() == comp_vec->NComps());
+   }
+   else
+   {
+      DBG_ASSERT(NComps_Cols() == 1);
+   }
+#endif
+   if( comp_vec )
+   {
+      if( NComps_Cols() != comp_vec->NComps() )
+      {
+         comp_vec = NULL;
+      }
+   }
+
+   for( Index irow = 0; irow < NComps_Rows(); irow++ )
+   {
+      for( Index jcol = 0; jcol < NComps_Cols(); jcol++ )
+      {
+         if( ConstComp(irow, jcol) )
+         {
+            SmartPtr<Vector> vec_j;
+            if( comp_vec )
+            {
+               vec_j = comp_vec->GetCompNonConst(jcol);
+            }
+            else
+            {
+               vec_j = &cols_norms;
+            }
+            DBG_ASSERT(IsValid(vec_j));
+            ConstComp(irow, jcol)->ComputeColA1(*vec_j, false);
+         }
+      }
+   }
+}
+
 void CompoundMatrix::PrintImpl(
    const Journalist&  jnlst,
    EJournalLevel      level,

src/LinAlg/IpCompoundMatrix.hpp

@@ -155,6 +155,16 @@ class IPOPTLIB_EXPORT CompoundMatrix: public Matrix
       bool    init
    ) const;
 
+   virtual void ComputeRowA1Impl(
+      Vector& rows_norms,
+      bool    init
+   ) const;
+
+   virtual void ComputeColA1Impl(
+      Vector& cols_norms,
+      bool    init
+   ) const;
+
    virtual void PrintImpl(
       const Journalist&  jnlst,
       EJournalLevel      level,

src/LinAlg/IpCompoundSymMatrix.cpp

@@ -237,6 +237,57 @@ void CompoundSymMatrix::ComputeRowAMaxImpl(
    }
 }
 
+void CompoundSymMatrix::ComputeRowA1Impl(
+   Vector& rows_norms,
+   bool    /*init*/
+) const
+{
+   if( !matrices_valid_ )
+   {
+      matrices_valid_ = MatricesValid();
+   }
+   DBG_ASSERT(matrices_valid_);
+
+   // The vector is assumed to be compound Vectors as well except if
+   // there is only one component
+   CompoundVector* comp_vec = dynamic_cast<CompoundVector*>(&rows_norms);
+
+   //  A few sanity checks
+   if( comp_vec )
+   {
+      DBG_ASSERT(NComps_Dim() == comp_vec->NComps());
+   }
+   else
+   {
+      DBG_ASSERT(NComps_Dim() == 1);
+   }
+
+   for( Index jcol = 0; jcol < NComps_Dim(); jcol++ )
+   {
+      for( Index irow = 0; irow < NComps_Dim(); irow++ )
+      {
+         SmartPtr<Vector> vec_i;
+         if( comp_vec )
+         {
+            vec_i = comp_vec->GetCompNonConst(irow);
+         }
+         else
+         {
+            vec_i = &rows_norms;
+         }
+         DBG_ASSERT(IsValid(vec_i));
+         if( jcol <= irow && ConstComp(irow, jcol) )
+         {
+            ConstComp(irow, jcol)->ComputeRowA1(*vec_i, false);
+         }
+         else if( jcol > irow && ConstComp(jcol, irow) )
+         {
+            ConstComp(jcol, irow)->ComputeRowA1(*vec_i, false);
+         }
+      }
+   }
+}
+
 void CompoundSymMatrix::PrintImpl(
    const Journalist&  jnlst,
    EJournalLevel      level,

src/LinAlg/IpCompoundSymMatrix.hpp

@@ -117,6 +117,11 @@ class IPOPTLIB_EXPORT CompoundSymMatrix: public SymMatrix
       bool    init
    ) const;
 
+   virtual void ComputeRowA1Impl(
+      Vector& rows_norms,
+      bool    init
+   ) const;
+
    virtual void PrintImpl(
       const Journalist&  jnlst,
       EJournalLevel      level,

src/LinAlg/IpDenseGenMatrix.cpp

@@ -415,6 +415,22 @@ void DenseGenMatrix::ComputeColAMaxImpl(
    }
 }
 
+void DenseGenMatrix::ComputeRowA1Impl(
+   Vector& /*rows_norms*/,
+   bool    /*init*/
+) const
+{
+   THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "DenseGenMatrix::ComputeRowA1Impl not implemented");
+}
+
+void DenseGenMatrix::ComputeColA1Impl(
+   Vector& /*cols_norms*/,
+   bool    /*init*/
+) const
+{
+   THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "DenseGenMatrix::ComputeColA1Impl not implemented");
+}
+
 bool DenseGenMatrix::HasValidNumbersImpl() const
 {
    DBG_ASSERT(initialized_);

src/LinAlg/IpDenseGenMatrix.hpp

@@ -217,6 +217,16 @@ class IPOPTLIB_EXPORT DenseGenMatrix: public Matrix
       bool    init
    ) const;
 
+   virtual void ComputeRowA1Impl(
+      Vector& rows_norms,
+      bool    init
+   ) const;
+
+   virtual void ComputeColA1Impl(
+      Vector& cols_norms,
+      bool    init
+   ) const;
+
    virtual void PrintImpl(
       const Journalist&  jnlst,
       EJournalLevel      level,

src/LinAlg/IpDenseSymMatrix.cpp

@@ -254,6 +254,15 @@ void DenseSymMatrix::ComputeRowAMaxImpl(
    }
 }
 
+void DenseSymMatrix::ComputeRowA1Impl(
+   Vector& /*rows_norms*/,
+   bool    /*init*/
+) const
+{
+   THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "DenseSymMatrix::ComputeRowA1Impl not implemented");
+}
+
+
 void DenseSymMatrix::PrintImpl(
    const Journalist&  jnlst,
    EJournalLevel      level,

src/LinAlg/IpDenseSymMatrix.hpp

@@ -133,6 +133,11 @@ class IPOPTLIB_EXPORT DenseSymMatrix: public SymMatrix
       bool    init
    ) const;
 
+   virtual void ComputeRowA1Impl(
+      Vector& rows_norms,
+      bool    init
+   ) const;
+
    virtual void PrintImpl(
       const Journalist&  jnlst,
       EJournalLevel      level,

src/LinAlg/IpDiagMatrix.cpp

@@ -71,6 +71,25 @@ void DiagMatrix::ComputeRowAMaxImpl(
    }
 }
 
+void DiagMatrix::ComputeRowA1Impl(
+   Vector& rows_norms,
+   bool    init
+) const
+{
+   DBG_ASSERT(IsValid(diag_));
+   if( init )
+   {
+      rows_norms.Copy(*diag_);
+      rows_norms.ElementWiseAbs();
+   }
+   else
+   {
+      SmartPtr<Vector> v = diag_->MakeNewCopy();
+      v->ElementWiseAbs();
+      rows_norms.Axpy(Number(1.), *v);
+   }
+}
+
 void DiagMatrix::PrintImpl(
    const Journalist&  jnlst,
    EJournalLevel      level,

src/LinAlg/IpDiagMatrix.hpp

@@ -63,6 +63,11 @@ class IPOPTLIB_EXPORT DiagMatrix: public SymMatrix
       bool    init
    ) const;
 
+   virtual void ComputeRowA1Impl(
+      Vector& rows_norms,
+      bool    init
+   ) const;
+
    virtual void PrintImpl(
       const Journalist&  jnlst,
       EJournalLevel      level,