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  spxgeometsc.hpp
26   	 * @brief Geometric mean row/column scaling.
27   	 */
28   	#include <assert.h>
29   	
30   	#include "soplex/spxout.h"
31   	#include "soplex/spxlpbase.h"
32   	#include "soplex/spxequilisc.h"
33   	
34   	namespace soplex
35   	{
36   	
37   	template <class R>
38   	static R computeScalingVec(
39   	   const SVSetBase<R>*   vecset, const std::vector<R>& coScaleval,
40   	   std::vector<R>&       scaleval, R epsilon)
41   	{
42   	   R pmax = 0.0;
43   	
44   	   assert(scaleval.size() == unsigned(vecset->num()));
45   	
46   	   for(int i = 0; i < vecset->num(); ++i)
47   	   {
48   	      const SVectorBase<R>& vec = (*vecset)[i];
49   	
50   	      R maxi = 0.0;
51   	      R mini = R(infinity);
52   	
53   	      for(int j = 0; j < vec.size(); ++j)
54   	      {
55   	         const R x = spxAbs(vec.value(j) * coScaleval[unsigned(vec.index(j))]);
56   	
57   	         if(!isZero(x, epsilon))
58   	         {
59   	            if(x > maxi)
60   	               maxi = x;
61   	
62   	            if(x < mini)
63   	               mini = x;
64   	         }
65   	      }
66   	
67   	      // empty rows/cols are possible
68   	      if(mini == R(infinity) || maxi == 0.0)
69   	      {
70   	         mini = 1.0;
71   	         maxi = 1.0;
72   	      }
73   	
74   	      assert(mini < R(infinity));
75   	      assert(maxi > 0.0);
76   	
77   	      scaleval[unsigned(i)] = 1.0 / spxSqrt(mini * maxi);
78   	
79   	      const R p = maxi / mini;
80   	
81   	      if(p > pmax)
82   	         pmax = p;
83   	   }
84   	
85   	   return pmax;
86   	}
87   	
88   	template <class R>
89   	SPxGeometSC<R>::SPxGeometSC(bool equilibrate, int maxIters, R minImpr, R goodEnough)
90   	   : SPxScaler<R>("Geometric")
91   	   , postequilibration(equilibrate)
92   	   , m_maxIterations(maxIters)
93   	   , m_minImprovement(minImpr)
94   	   , m_goodEnoughRatio(goodEnough)
95   	{
96   	   assert(maxIters > 0);
97   	   assert(minImpr > 0.0 && minImpr <= 1.0);
98   	   assert(goodEnough >= 0.0);
99   	}
100  	
101  	template <class R>
102  	SPxGeometSC<R>::SPxGeometSC(const SPxGeometSC<R>& old)
103  	   : SPxScaler<R>(old)
104  	   , postequilibration(old.postequilibration)
105  	   , m_maxIterations(old.m_maxIterations)
106  	   , m_minImprovement(old.m_minImprovement)
107  	   , m_goodEnoughRatio(old.m_goodEnoughRatio)
108  	{
109  	   assert(m_maxIterations > 0);
110  	   assert(m_minImprovement > 0.0 && m_minImprovement <= 1.0);
111  	   assert(m_goodEnoughRatio >= 0.0);
112  	}
113  	
114  	template <class R>
115  	SPxGeometSC<R>& SPxGeometSC<R>::operator=(const SPxGeometSC<R>& rhs)
116  	{
117  	   if(this != &rhs)
118  	   {
119  	      SPxScaler<R>::operator=(rhs);
120  	   }
121  	
122  	   return *this;
123  	}
124  	
125  	template <class R>
126  	void SPxGeometSC<R>::scale(SPxLPBase<R>& lp, bool persistent)
127  	{
128  	
129  	   SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "Geometric scaling LP" <<
130  	                 (persistent ? " (persistent)" : "") << (postequilibration ? " with post-equilibration" : "") <<
131  	                 std::endl;)
132  	
133  	   this->setup(lp);
134  	
135  	   /* We want to do that direction first, with the lower ratio.
136  	    * See SPxEquiliSC<R>::scale() for a reasoning.
137  	    */
138  	   const R colratio = this->maxColRatio(lp);
139  	   const R rowratio = this->maxRowRatio(lp);
140  	   R epsilon = this->tolerances()->epsilon();
141  	
142  	   const bool colFirst = colratio < rowratio;
143  	
144  	   R p0start;
145  	   R p1start;
146  	
147  	   if(colFirst)
148  	   {
149  	      p0start = colratio;
150  	      p1start = rowratio;
151  	   }
152  	   else
153  	   {
154  	      p0start = rowratio;
155  	      p1start = colratio;
156  	   }
157  	
158  	   SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "before scaling:"
159  	                 << " min= " << lp.minAbsNzo()
160  	                 << " max= " << lp.maxAbsNzo()
161  	                 << " col-ratio= " << colratio
162  	                 << " row-ratio= " << rowratio
163  	                 << std::endl;)
164  	
165  	   // perform geometric scaling only if maximum ratio is above threshold
166  	   bool geoscale = p1start > m_goodEnoughRatio;
167  	
168  	   if(!geoscale)
169  	   {
170  	      SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "No geometric scaling done, ratio good enough" <<
171  	                    std::endl;)
172  	
173  	      if(!postequilibration)
174  	      {
175  	         lp.setScalingInfo(true);
176  	         return;
177  	      }
178  	
179  	      SPX_MSG_INFO2((*this->spxout),
180  	                    (*this->spxout) << " ... but will still perform equilibrium scaling" <<
181  	                    std::endl;)
182  	   }
183  	
184  	   std::vector<R> rowscale(unsigned(lp.nRows()), 1.0);
185  	   std::vector<R> colscale(unsigned(lp.nCols()), 1.0);
186  	
187  	   R p0 = 0.0;
188  	   R p1 = 0.0;
189  	
190  	   if(geoscale)
191  	   {
192  	      R p0prev = p0start;
193  	      R p1prev = p1start;
194  	
195  	      // we make at most maxIterations.
196  	      for(int count = 0; count < m_maxIterations; count++)
197  	      {
198  	         if(colFirst)
199  	         {
200  	            p0 = computeScalingVec(lp.colSet(), rowscale, colscale, epsilon);
201  	            p1 = computeScalingVec(lp.rowSet(), colscale, rowscale, epsilon);
202  	         }
203  	         else
204  	         {
205  	            p0 = computeScalingVec(lp.rowSet(), colscale, rowscale, epsilon);
206  	            p1 = computeScalingVec(lp.colSet(), rowscale, colscale, epsilon);
207  	         }
208  	
209  	         SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "Geometric scaling round " << count
210  	                       << " col-ratio= " << (colFirst ? p0 : p1)
211  	                       << " row-ratio= " << (colFirst ? p1 : p0)
212  	                       << std::endl;)
213  	
214  	         if(p0 > m_minImprovement * p0prev && p1 > m_minImprovement * p1prev)
215  	            break;
216  	
217  	         p0prev = p0;
218  	         p1prev = p1;
219  	      }
220  	
221  	      // perform geometric scaling only if there is enough (default 15%) improvement.
222  	      geoscale = (p0 <= m_minImprovement * p0start || p1 <= m_minImprovement * p1start);
223  	   }
224  	
225  	   if(!geoscale && !postequilibration)
226  	   {
227  	      SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "No geometric scaling done." << std::endl;)
228  	      lp.setScalingInfo(true);
229  	   }
230  	   else
231  	   {
232  	      DataArray<int>& colscaleExp = *this->m_activeColscaleExp;
233  	      DataArray<int>& rowscaleExp = *this->m_activeRowscaleExp;
234  	
235  	      if(postequilibration)
236  	      {
237  	         if(!geoscale)
238  	         {
239  	            std::fill(rowscale.begin(), rowscale.end(), 1.0);
240  	            std::fill(colscale.begin(), colscale.end(), 1.0);
241  	         }
242  	
243  	         SPxEquiliSC<R>::computePostequiExpVecs(lp, rowscale, colscale, rowscaleExp, colscaleExp, epsilon);
244  	      }
245  	      else
246  	      {
247  	         this->computeExpVec(colscale, colscaleExp);
248  	         this->computeExpVec(rowscale, rowscaleExp);
249  	      }
250  	
251  	      this->applyScaling(lp);
252  	
253  	      SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "Row scaling min= " << this->minAbsRowscale()
254  	                    << " max= " << this->maxAbsRowscale()
255  	                    << std::endl
256  	                    << "IGEOSC06 Col scaling min= " << this->minAbsColscale()
257  	                    << " max= " << this->maxAbsColscale()
258  	                    << std::endl;)
259  	
260  	      SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "after scaling: "
261  	                    << " min= " << lp.minAbsNzo(false)
262  	                    << " max= " << lp.maxAbsNzo(false)
263  	                    << " col-ratio= " << this->maxColRatio(lp)
264  	                    << " row-ratio= " << this->maxRowRatio(lp)
265  	                    << std::endl;)
266  	   }
267  	}
268  	
269  	
270  	} // namespace soplex
271