.. _program_listing_file__tmp_ws_src_proxsuite_include_proxsuite_proxqp_dense_preconditioner_ruiz.hpp:

Program Listing for File ruiz.hpp
=================================

|exhale_lsh| :ref:`Return to documentation for file <file__tmp_ws_src_proxsuite_include_proxsuite_proxqp_dense_preconditioner_ruiz.hpp>` (``/tmp/ws/src/proxsuite/include/proxsuite/proxqp/dense/preconditioner/ruiz.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   //
   // Copyright (c) 2022 INRIA
   //
   #ifndef PROXSUITE_PROXQP_DENSE_PRECOND_RUIZ_HPP
   #define PROXSUITE_PROXQP_DENSE_PRECOND_RUIZ_HPP
   
   #include "proxsuite/proxqp/dense/views.hpp"
   #include "proxsuite/proxqp/dense/fwd.hpp"
   #include <proxsuite/linalg/dense/core.hpp>
   #include <proxsuite/proxqp/settings.hpp>
   #include <ostream>
   
   #include <Eigen/Core>
   
   namespace proxsuite {
   namespace proxqp {
   enum struct Symmetry
   {
     general,
     lower,
     upper,
   };
   namespace dense {
   namespace detail {
   
   template<typename T>
   auto
   ruiz_scale_qp_in_place( //
     VectorViewMut<T> delta_,
     std::ostream* logger_ptr,
     QpViewBoxMut<T> qp,
     T epsilon,
     isize max_iter,
     Symmetry sym,
     ProblemType problem_type,
     proxsuite::linalg::veg::dynstack::DynStackMut stack) -> T
   {
   
     T c(1);
     auto S = delta_.to_eigen();
     auto H = qp.H.to_eigen();
     auto g = qp.g.to_eigen();
     auto A = qp.A.to_eigen();
     auto b = qp.b.to_eigen();
     auto C = qp.C.to_eigen();
     auto u = qp.u.to_eigen();
     auto l = qp.l.to_eigen();
   
     static constexpr T machine_eps = std::numeric_limits<T>::epsilon();
     /*
      * compute equilibration parameters and scale in place the qp following
      * algorithm
      *
      * modified: removed g in gamma computation
      */
   
     isize n = qp.H.rows;
     isize n_eq = qp.A.rows;
     isize n_in = qp.C.rows;
   
     T gamma = T(1);
   
     LDLT_TEMP_VEC(T, delta, n + n_eq + n_in, stack);
   
     i64 iter = 1;
     while (infty_norm((1 - delta.array()).matrix()) > epsilon) {
       if (logger_ptr != nullptr) {
         *logger_ptr                                   //
           << "j : "                                   //
           << iter                                     //
           << " ; error : "                            //
           << infty_norm((1 - delta.array()).matrix()) //
           << "\n\n";
       }
       if (iter == max_iter) {
         break;
       } else {
         ++iter;
       }
       // normalization vector
       {
         switch (problem_type) {
           case ProblemType::LP:
             for (isize k = 0; k < n; ++k) {
               T aux = sqrt(std::max({
                 n_eq > 0 ? infty_norm(A.col(k)) : T(0),
                 n_in > 0 ? infty_norm(C.col(k)) : T(0),
               }));
               if (aux == T(0)) {
                 delta(k) = T(1);
               } else {
                 delta(k) = T(1) / (aux + machine_eps);
               }
             }
             break;
           case ProblemType::QP:
             for (isize k = 0; k < n; ++k) {
               switch (sym) {
                 case Symmetry::upper: { // upper triangular part
                   T aux = sqrt(std::max({
                     infty_norm(H.col(k).head(k)),
                     infty_norm(H.row(k).tail(n - k)),
                     n_eq > 0 ? infty_norm(A.col(k)) : T(0),
                     n_in > 0 ? infty_norm(C.col(k)) : T(0),
                   }));
                   if (aux == T(0)) {
                     delta(k) = T(1);
                   } else {
                     delta(k) = T(1) / (aux + machine_eps);
                   }
                   break;
                 }
                 case Symmetry::lower: { // lower triangular part
   
                   T aux = sqrt(std::max({
                     infty_norm(H.col(k).head(k)),
                     infty_norm(H.col(k).tail(n - k)),
                     n_eq > 0 ? infty_norm(A.col(k)) : T(0),
                     n_in > 0 ? infty_norm(C.col(k)) : T(0),
                   }));
                   if (aux == T(0)) {
                     delta(k) = T(1);
                   } else {
                     delta(k) = T(1) / (aux + machine_eps);
                   }
                   break;
                 }
                 case Symmetry::general: {
   
                   T aux = sqrt(std::max({
                     infty_norm(H.col(k)),
                     n_eq > 0 ? infty_norm(A.col(k)) : T(0),
                     n_in > 0 ? infty_norm(C.col(k)) : T(0),
                   }));
                   if (aux == T(0)) {
                     delta(k) = T(1);
                   } else {
                     delta(k) = T(1) / (aux + machine_eps);
                   }
                   break;
                 }
               }
             }
             break;
         }
         for (isize k = 0; k < n_eq; ++k) {
           T aux = sqrt(infty_norm(A.row(k)));
           if (aux == T(0)) {
             delta(n + k) = T(1);
           } else {
             delta(n + k) = T(1) / (aux + machine_eps);
           }
         }
         for (isize k = 0; k < n_in; ++k) {
           T aux = sqrt(infty_norm(C.row(k)));
           if (aux == T(0)) {
             delta(k + n + n_eq) = T(1);
           } else {
             delta(k + n + n_eq) = T(1) / (aux + machine_eps);
           }
         }
       }
       {
   
         // normalize A and C
         A = delta.segment(n, n_eq).asDiagonal() * A * delta.head(n).asDiagonal();
         C = delta.tail(n_in).asDiagonal() * C * delta.head(n).asDiagonal();
         // normalize vectors
         g.array() *= delta.head(n).array();
         b.array() *= delta.middleRows(n, n_eq).array();
         u.array() *= delta.tail(n_in).array();
         l.array() *= delta.tail(n_in).array();
   
         // normalize H
         switch (problem_type) {
           case ProblemType::LP:
             break;
           case ProblemType::QP:
             switch (sym) {
               case Symmetry::upper: {
                 // upper triangular part
                 for (isize j = 0; j < n; ++j) {
                   H.col(j).head(j + 1) *= delta(j);
                 }
                 // normalisation des lignes
                 for (isize i = 0; i < n; ++i) {
                   H.row(i).tail(n - i) *= delta(i);
                 }
                 break;
               }
               case Symmetry::lower: {
                 // lower triangular part
                 for (isize j = 0; j < n; ++j) {
                   H.col(j).tail(n - j) *= delta(j);
                 }
                 // normalisation des lignes
                 for (isize i = 0; i < n; ++i) {
                   H.row(i).head(i + 1) *= delta(i);
                 }
                 break;
               }
               case Symmetry::general: {
                 // all matrix
                 H = delta.head(n).asDiagonal() * H * delta.head(n).asDiagonal();
                 break;
               }
               default:
                 break;
             }
   
             // additional normalization for the cost function
             switch (sym) {
               case Symmetry::upper: {
                 // upper triangular part
                 T tmp = T(0);
                 for (isize j = 0; j < n; ++j) {
                   tmp += proxqp::dense::infty_norm(H.row(j).tail(n - j));
                 }
                 gamma = 1 / std::max(tmp / T(n), T(1));
                 break;
               }
               case Symmetry::lower: {
                 // lower triangular part
                 T tmp = T(0);
                 for (isize j = 0; j < n; ++j) {
                   tmp += proxqp::dense::infty_norm(H.col(j).tail(n - j));
                 }
                 gamma = 1 / std::max(tmp / T(n), T(1));
                 break;
               }
               case Symmetry::general: {
                 // all matrix
                 gamma =
                   1 / std::max(
                         T(1),
                         (H.colwise().template lpNorm<Eigen::Infinity>()).mean());
                 break;
               }
               default:
                 break;
             }
             H *= gamma;
             break;
         }
         g *= gamma;
   
         S.array() *= delta.array(); // coefficientwise product
         c *= gamma;
       }
     }
     return c;
   }
   } // namespace detail
   
   namespace preconditioner {
   
   template<typename T>
   struct RuizEquilibration
   {
     Vec<T> delta;
     T c;
     isize dim;
     T epsilon;
     i64 max_iter;
     Symmetry sym;
   
     std::ostream* logger_ptr = nullptr;
     explicit RuizEquilibration(isize dim_,
                                isize n_eq_in,
                                T epsilon_ = T(1e-3),
                                i64 max_iter_ = 10,
                                Symmetry sym_ = Symmetry::general,
                                std::ostream* logger = nullptr)
       : delta(Vec<T>::Ones(dim_ + n_eq_in))
       , c(1)
       , dim(dim_)
       , epsilon(epsilon_)
       , max_iter(max_iter_)
       , sym(sym_)
       , logger_ptr(logger)
     {
     }
     void print() const
     {
       // CHANGE: endl to newline
       *logger_ptr << " delta : " << delta << "\n\n";
       *logger_ptr << " c : " << c << "\n\n";
     }
     static auto scale_qp_in_place_req(proxsuite::linalg::veg::Tag<T> tag,
                                       isize n,
                                       isize n_eq,
                                       isize n_in)
       -> proxsuite::linalg::veg::dynstack::StackReq
     {
       return proxsuite::linalg::dense::temp_vec_req(tag, n + n_eq + n_in);
     }
   
     // H_new = c * head @ H @ head
     // A_new = tail @ A @ head
     // g_new = c * head @ g
     // b_new = tail @ b
     void scale_qp_in_place(QpViewBoxMut<T> qp,
                            bool execute_preconditioner,
                            const isize max_iter,
                            const T epsilon,
                            const ProblemType problem_type,
                            proxsuite::linalg::veg::dynstack::DynStackMut stack)
     {
       if (execute_preconditioner) {
         delta.setOnes();
         c = detail::ruiz_scale_qp_in_place({ proxqp::from_eigen, delta },
                                            logger_ptr,
                                            qp,
                                            epsilon,
                                            max_iter,
                                            sym,
                                            problem_type,
                                            stack);
       } else {
   
         auto H = qp.H.to_eigen();
         auto g = qp.g.to_eigen();
         auto A = qp.A.to_eigen();
         auto b = qp.b.to_eigen();
         auto C = qp.C.to_eigen();
         auto u = qp.u.to_eigen();
         auto l = qp.l.to_eigen();
         isize n = qp.H.rows;
         isize n_eq = qp.A.rows;
         isize n_in = qp.C.rows;
   
         // normalize A and C
         A = delta.segment(n, n_eq).asDiagonal() * A * delta.head(n).asDiagonal();
         C = delta.tail(n_in).asDiagonal() * C * delta.head(n).asDiagonal();
   
         // normalize H
         switch (sym) {
           case Symmetry::upper: {
             // upper triangular part
             for (isize j = 0; j < n; ++j) {
               H.col(j).head(j + 1) *= delta(j);
             }
             // normalisation des lignes
             for (isize i = 0; i < n; ++i) {
               H.row(i).tail(n - i) *= delta(i);
             }
             break;
           }
           case Symmetry::lower: {
             // lower triangular part
             for (isize j = 0; j < n; ++j) {
               H.col(j).tail(n - j) *= delta(j);
             }
             // normalisation des lignes
             for (isize i = 0; i < n; ++i) {
               H.row(i).head(i + 1) *= delta(i);
             }
             break;
           }
           case Symmetry::general: {
             // all matrix
             H = delta.head(n).asDiagonal() * H * delta.head(n).asDiagonal();
             break;
           }
           default:
             break;
         }
   
         // normalize vectors
         g.array() *= delta.head(n).array();
         b.array() *= delta.segment(n, n_eq).array();
         l.array() *= delta.tail(n_in).array();
         u.array() *= delta.tail(n_in).array();
   
         g *= c;
         H *= c;
       }
     }
     void scale_qp(const QpViewBox<T> qp,
                   QpViewBoxMut<T> scaled_qp,
                   bool execute_preconditioner,
                   VectorViewMut<T> tmp_delta_preallocated) const
     {
   
       /*
        * scaled_qp is scaled, whereas first qp is not
        * the procedure computes as well equilibration parameters using default
        * parameters
        */
   
       scaled_qp.H.to_eigen() = qp.H.to_eigen();
       scaled_qp.A.to_eigen() = qp.A.to_eigen();
       scaled_qp.C.to_eigen() = qp.C.to_eigen();
       scaled_qp.g.to_eigen() = qp.g.to_eigen();
       scaled_qp.b.to_eigen() = qp.b.to_eigen();
       scaled_qp.d.to_eigen() = qp.d.to_eigen();
   
       scale_qp_in_place(scaled_qp,
                         execute_preconditioner,
                         tmp_delta_preallocated,
                         epsilon,
                         max_iter);
     }
     // modifies variables in place
     void scale_primal_in_place(VectorViewMut<T> primal) const
     {
       primal.to_eigen().array() /= delta.array().head(dim);
     }
     void scale_dual_in_place(VectorViewMut<T> dual) const
     {
       dual.to_eigen().array() = dual.as_const().to_eigen().array() /
                                 delta.tail(delta.size() - dim).array() * c;
     }
     void scale_dual_in_place_eq(VectorViewMut<T> dual) const
     {
       dual.to_eigen().array() =
         dual.as_const().to_eigen().array() /
         delta.middleRows(dim, dual.to_eigen().size()).array() * c;
     }
     void scale_dual_in_place_in(VectorViewMut<T> dual) const
     {
       dual.to_eigen().array() = dual.as_const().to_eigen().array() /
                                 delta.tail(dual.to_eigen().size()).array() * c;
     }
     void unscale_primal_in_place(VectorViewMut<T> primal) const
     {
       primal.to_eigen().array() *= delta.array().head(dim);
     }
     void unscale_dual_in_place(VectorViewMut<T> dual) const
     {
       dual.to_eigen().array() = dual.as_const().to_eigen().array() *
                                 delta.tail(delta.size() - dim).array() / c;
     }
     void unscale_dual_in_place_eq(VectorViewMut<T> dual) const
     {
       dual.to_eigen().array() =
         dual.as_const().to_eigen().array() *
         delta.middleRows(dim, dual.to_eigen().size()).array() / c;
     }
     void unscale_dual_in_place_in(VectorViewMut<T> dual) const
     {
       dual.to_eigen().array() = dual.as_const().to_eigen().array() *
                                 delta.tail(dual.to_eigen().size()).array() / c;
     }
     // modifies residuals in place
     void scale_primal_residual_in_place(VectorViewMut<T> primal) const
     {
       primal.to_eigen().array() *= delta.tail(delta.size() - dim).array();
     }
   
     void scale_primal_residual_in_place_eq(VectorViewMut<T> primal_eq) const
     {
       primal_eq.to_eigen().array() *=
         delta.middleRows(dim, primal_eq.to_eigen().size()).array();
     }
     void scale_primal_residual_in_place_in(VectorViewMut<T> primal_in) const
     {
       primal_in.to_eigen().array() *=
         delta.tail(primal_in.to_eigen().size()).array();
     }
     void scale_dual_residual_in_place(VectorViewMut<T> dual) const
     {
       dual.to_eigen().array() *= delta.head(dim).array() * c;
     }
     void unscale_primal_residual_in_place(VectorViewMut<T> primal) const
     {
       primal.to_eigen().array() /= delta.tail(delta.size() - dim).array();
     }
     void unscale_primal_residual_in_place_eq(VectorViewMut<T> primal_eq) const
     {
       primal_eq.to_eigen().array() /=
         delta.middleRows(dim, primal_eq.to_eigen().size()).array();
     }
     void unscale_primal_residual_in_place_in(VectorViewMut<T> primal_in) const
     {
       primal_in.to_eigen().array() /=
         delta.tail(primal_in.to_eigen().size()).array();
     }
     void unscale_dual_residual_in_place(VectorViewMut<T> dual) const
     {
       dual.to_eigen().array() /= delta.head(dim).array() * c;
     }
   };
   
   } // namespace preconditioner
   } // namespace dense
   } // namespace proxqp
   } // namespace proxsuite
   
   #endif /* end of include guard PROXSUITE_PROXQP_DENSE_PRECOND_RUIZ_HPP */