1    	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2    	/*                                                                           */
3    	/*                  This file is part of the class library                   */
4    	/*       SoPlex --- the Sequential object-oriented simPlex.                  */
5    	/*                                                                           */
6    	/*  Copyright (c) 1996-2023 Zuse Institute Berlin (ZIB)                      */
7    	/*                                                                           */
8    	/*  Licensed under the Apache License, Version 2.0 (the "License");          */
9    	/*  you may not use this file except in compliance with the License.         */
10   	/*  You may obtain a copy of the License at                                  */
11   	/*                                                                           */
12   	/*      http://www.apache.org/licenses/LICENSE-2.0                           */
13   	/*                                                                           */
14   	/*  Unless required by applicable law or agreed to in writing, software      */
15   	/*  distributed under the License is distributed on an "AS IS" BASIS,        */
16   	/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17   	/*  See the License for the specific language governing permissions and      */
18   	/*  limitations under the License.                                           */
19   	/*                                                                           */
20   	/*  You should have received a copy of the Apache-2.0 license                */
21   	/*  along with SoPlex; see the file LICENSE. If not email to soplex@zib.de.  */
22   	/*                                                                           */
23   	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24   	
25   	/**@file  spxequilisc.hpp
26   	 * @brief Equilibrium row/column scaling.
27   	 */
28   	#include <assert.h>
29   	
30   	#include "soplex/spxout.h"
31   	#include "soplex/spxlpbase.h"
32   	#include "soplex/spxlp.h"
33   	#include "soplex.h"
34   	
35   	namespace soplex
36   	{
37   	static inline const char* makename(bool doBoth)
38   	{
39   	   return doBoth ? "bi-Equilibrium" : "uni-Equilibrium";
40   	}
41   	
42   	/// maximum ratio between absolute biggest and smallest element in any scaled row/column.
43   	template <class R>
44   	static R maxPrescaledRatio(const SPxLPBase<R>& lp, const std::vector<R>& coScaleval, bool rowRatio)
45   	{
46   	   R pmax = 0.0;
47   	   const int n = rowRatio ? lp.nRows() : lp.nCols();
48   	
49   	   for(int i = 0; i < n; ++i)
50   	   {
51   	      const SVectorBase<R>& vec = rowRatio ? lp.rowVector(i) : lp.colVector(i);
52   	      R mini = R(infinity);
53   	      R maxi = 0.0;
54   	
55   	      for(int j = 0; j < vec.size(); ++j)
56   	      {
57   	         assert(vec.index(j) >= 0);
58   	         const R x = spxAbs(vec.value(j)) * coScaleval[unsigned(vec.index(j))];
59   	
60   	         if(isZero(x, lp.tolerances()->epsilon()))
61   	            continue;
62   	
63   	         if(x < mini)
64   	            mini = x;
65   	
66   	         if(x > maxi)
67   	            maxi = x;
68   	      }
69   	
70   	      if(mini == R(infinity))
71   	         continue;
72   	
73   	      const R p = maxi / mini;
74   	
75   	      if(p > pmax)
76   	         pmax = p;
77   	   }
78   	
79   	   return pmax;
80   	}
81   	
82   	template <class R>
83   	void SPxEquiliSC<R>::computeEquiExpVec(const SVSetBase<R>* vecset, const DataArray<int>& coScaleExp,
84   	                                       DataArray<int>& scaleExp, R epsilon)
85   	{
86   	   assert(vecset != nullptr);
87   	
88   	   for(int i = 0; i < vecset->num(); ++i)
89   	   {
90   	      const SVectorBase<R>& vec = (*vecset)[i];
91   	
92   	      R maxi = 0.0;
93   	
94   	      for(int j = 0; j < vec.size(); ++j)
95   	      {
96   	         const R x = spxAbs(spxLdexp(vec.value(j), coScaleExp[vec.index(j)]));
97   	
98   	         if(GT(x, maxi, epsilon))
99   	            maxi = x;
100  	      }
101  	
102  	      // empty rows/cols are possible
103  	      if(maxi == 0.0)
104  	         maxi = 1.0;
105  	
106  	      assert(maxi > 0.0);
107  	
108  	      spxFrexp(Real(1.0 / maxi), &(scaleExp[i]));
109  	
110  	      scaleExp[i] -= 1;
111  	   }
112  	}
113  	
114  	template <class R>
115  	void SPxEquiliSC<R>::computeEquiExpVec(const SVSetBase<R>* vecset, const std::vector<R>& coScaleVal,
116  	                                       DataArray<int>& scaleExp, R epsilon)
117  	{
118  	   assert(vecset != nullptr);
119  	
120  	   for(int i = 0; i < vecset->num(); ++i)
121  	   {
122  	      const SVectorBase<R>& vec = (*vecset)[i];
123  	
124  	      R maxi = 0.0;
125  	
126  	      for(int j = 0; j < vec.size(); ++j)
127  	      {
128  	         assert(vec.index(j) >= 0);
129  	         const R x = spxAbs(vec.value(j) * coScaleVal[unsigned(vec.index(j))]);
130  	
131  	         if(GT(x, maxi, epsilon))
132  	            maxi = x;
133  	      }
134  	
135  	      // empty rows/cols are possible
136  	      if(maxi == 0.0)
137  	         maxi = 1.0;
138  	
139  	      assert(maxi > 0.0);
140  	
141  	      spxFrexp(Real(1.0 / maxi), &(scaleExp[i]));
142  	
143  	      scaleExp[i] -= 1;
144  	   }
145  	}
146  	
147  	template <class R>
148  	void SPxEquiliSC<R>::computePostequiExpVecs(const SPxLPBase<R>& lp,
149  	      const std::vector<R>& preRowscale, const std::vector<R>& preColscale,
150  	      DataArray<int>& rowscaleExp, DataArray<int>& colscaleExp, R epsilon)
151  	{
152  	   const R colratio = maxPrescaledRatio(lp, preRowscale, false);
153  	   const R rowratio = maxPrescaledRatio(lp, preColscale, true);
154  	
155  	   const bool colFirst = colratio < rowratio;
156  	
157  	   // see SPxEquiliSC<R>::scale for reason behind this branch
158  	   if(colFirst)
159  	   {
160  	      computeEquiExpVec(lp.colSet(), preRowscale, colscaleExp, epsilon);
161  	      computeEquiExpVec(lp.rowSet(), colscaleExp, rowscaleExp, epsilon);
162  	   }
163  	   else
164  	   {
165  	      computeEquiExpVec(lp.rowSet(), preColscale, rowscaleExp, epsilon);
166  	      computeEquiExpVec(lp.colSet(), rowscaleExp, colscaleExp, epsilon);
167  	   }
168  	}
169  	
170  	template <class R>
171  	SPxEquiliSC<R>::SPxEquiliSC(bool doBoth)
172  	   : SPxScaler<R>(makename(doBoth), false, doBoth)
173  	{}
174  	
175  	template <class R>
176  	SPxEquiliSC<R>::SPxEquiliSC(const SPxEquiliSC<R>& old)
177  	   : SPxScaler<R>(old)
178  	{}
179  	
180  	template <class R>
181  	SPxEquiliSC<R>& SPxEquiliSC<R>::operator=(const SPxEquiliSC<R>& rhs)
182  	{
183  	   if(this != &rhs)
184  	   {
185  	      SPxScaler<R>::operator=(rhs);
186  	   }
187  	
188  	   return *this;
189  	}
190  	
191  	template <class R>
192  	void SPxEquiliSC<R>::scale(SPxLPBase<R>& lp, bool persistent)
193  	{
194  	
195  	   SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "Equilibrium scaling LP" <<
196  	                 (persistent ? " (persistent)" : "") << std::endl;)
197  	
198  	   this->setup(lp);
199  	
200  	   /* We want to do the direction first, which has a lower maximal ratio,
201  	    * since the lowest value in the scaled matrix is bounded from below by
202  	    * the inverse of the maximum ratio of the direction that is done first
203  	    * Example:
204  	    *
205  	    *                     Rowratio
206  	    *            0.1  1   10
207  	    *            10   1   10
208  	    *
209  	    * Colratio   100  1
210  	    *
211  	    * Row first =>         Col next =>
212  	    *            0.1  1          0.1  1
213  	    *            1    0.1        1    0.1
214  	    *
215  	    * Col first =>         Row next =>
216  	    *            0.01 1          0.01 1
217  	    *            1    1          1    1
218  	    *
219  	    */
220  	   R colratio = this->maxColRatio(lp);
221  	   R rowratio = this->maxRowRatio(lp);
222  	   R epsilon = this->tolerances()->epsilon();
223  	
224  	   bool colFirst = colratio < rowratio;
225  	
226  	   SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "before scaling:"
227  	                 << " min= " << lp.minAbsNzo()
228  	                 << " max= " << lp.maxAbsNzo()
229  	                 << " col-ratio= " << colratio
230  	                 << " row-ratio= " << rowratio
231  	                 << std::endl;)
232  	
233  	   if(colFirst)
234  	   {
235  	      computeEquiExpVec(lp.colSet(), *this->m_activeRowscaleExp, *this->m_activeColscaleExp, epsilon);
236  	
237  	      if(this->m_doBoth)
238  	         computeEquiExpVec(lp.rowSet(), *this->m_activeColscaleExp, *this->m_activeRowscaleExp, epsilon);
239  	   }
240  	   else
241  	   {
242  	      computeEquiExpVec(lp.rowSet(), *this->m_activeColscaleExp, *this->m_activeRowscaleExp, epsilon);
243  	
244  	      if(this->m_doBoth)
245  	         computeEquiExpVec(lp.colSet(), *this->m_activeRowscaleExp, *this->m_activeColscaleExp, epsilon);
246  	   }
247  	
248  	   /* scale */
249  	   this->applyScaling(lp);
250  	
251  	   SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "Row scaling min= " << this->minAbsRowscale()
252  	                 << " max= " << this->maxAbsRowscale()
253  	                 << std::endl
254  	                 << "Col scaling min= " << this->minAbsColscale()
255  	                 << " max= " << this->maxAbsColscale()
256  	                 << std::endl;)
257  	
258  	   SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "after scaling: "
259  	                 << " min= " << lp.minAbsNzo(false)
260  	                 << " max= " << lp.maxAbsNzo(false)
261  	                 << " col-ratio= " << this->maxColRatio(lp)
262  	                 << " row-ratio= " << this->maxRowRatio(lp)
263  	                 << std::endl;)
264  	
265  	}
266  	
267  	} // namespace soplex
268