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  spxscaler.hpp
26   	 * @brief LP scaling base class.
27   	 */
28   	
29   	#include <cmath>
30   	
31   	#include <iostream>
32   	#include <assert.h>
33   	#include "soplex/dsvector.h"
34   	#include "soplex/lprowsetbase.h"
35   	#include "soplex/lpcolsetbase.h"
36   	#include <limits>
37   	
38   	namespace soplex
39   	{
40   	
41   	template <class R>
42   	std::ostream& operator<<(std::ostream& s, const SPxScaler<R>& sc)
43   	{
44   	   const DataArray < int >& colscaleExp = *(sc.m_activeColscaleExp);
45   	   DataArray < int > rowccaleExp = *(sc.m_activeRowscaleExp);
46   	
47   	   s << sc.getName() << " scaler:" << std::endl;
48   	   s << "colscale = [ ";
49   	
50   	   for(int ci = 0; ci < colscaleExp.size(); ++ci)
51   	      s << colscaleExp[ci] << " ";
52   	
53   	   s << "]" << std::endl;
54   	
55   	   s << "rowscale = [ ";
56   	
57   	   for(int ri = 0; ri < rowccaleExp.size(); ++ri)
58   	      s << rowccaleExp[ri] << " ";
59   	
60   	   s << "]" << std::endl;
61   	
62   	   return s;
63   	}
64   	
65   	
66   	template <class R>
67   	SPxScaler<R>::SPxScaler(
68   	   const char* name,
69   	   bool        colFirst,
70   	   bool        doBoth,
71   	   SPxOut*     outstream)
72   	   : m_name(name)
73   	   , m_activeColscaleExp(0)
74   	   , m_activeRowscaleExp(0)
75   	   , m_colFirst(colFirst)
76   	   , m_doBoth(doBoth)
77   	   , spxout(outstream)
78   	{
79   	}
80   	
81   	template <class R>
82   	SPxScaler<R>::SPxScaler(const SPxScaler<R>& old)
83   	   : m_name(old.m_name)
84   	   , m_activeColscaleExp(old.m_activeColscaleExp)
85   	   , m_activeRowscaleExp(old.m_activeRowscaleExp)
86   	   , m_colFirst(old.m_colFirst)
87   	   , m_doBoth(old.m_doBoth)
88   	   , spxout(old.spxout)
89   	{
90   	}
91   	
92   	template <class R>
93   	SPxScaler<R>::~SPxScaler()
94   	{
95   	   m_name = 0;
96   	}
97   	
98   	template <class R>
99   	SPxScaler<R>& SPxScaler<R>::operator=(const SPxScaler<R>& rhs)
100  	{
101  	   if(this != &rhs)
102  	   {
103  	      m_name     = rhs.m_name;
104  	      m_activeColscaleExp = rhs.m_activeColscaleExp;
105  	      m_activeRowscaleExp = rhs.m_activeRowscaleExp;
106  	      m_colFirst = rhs.m_colFirst;
107  	      m_doBoth   = rhs.m_doBoth;
108  	      spxout     = rhs.spxout;
109  	
110  	      assert(SPxScaler<R>::isConsistent());
111  	   }
112  	
113  	   return *this;
114  	}
115  	
116  	template <class R>
117  	const char* SPxScaler<R>::getName() const
118  	{
119  	
120  	   return m_name;
121  	}
122  	
123  	template <class R>
124  	void SPxScaler<R>::setOrder(bool colFirst)
125  	{
126  	
127  	   m_colFirst = colFirst;
128  	}
129  	
130  	template <class R>
131  	void SPxScaler<R>::setBoth(bool both)
132  	{
133  	
134  	   m_doBoth = both;
135  	}
136  	
137  	template <class R>
138  	void SPxScaler<R>::setRealParam(R param, const char* name)
139  	{}
140  	
141  	template <class R>
142  	void SPxScaler<R>::setIntParam(int param, const char* name)
143  	{}
144  	
145  	template <class R>
146  	void SPxScaler<R>::setup(SPxLPBase<R>& lp)
147  	{
148  	   assert(lp.isConsistent());
149  	   m_activeColscaleExp = &lp.LPColSetBase<R>::scaleExp;
150  	   m_activeRowscaleExp = &lp.LPRowSetBase<R>::scaleExp;
151  	   m_activeColscaleExp->reSize(lp.nCols());
152  	   m_activeRowscaleExp->reSize(lp.nRows());
153  	
154  	   for(int i = 0; i < lp.nCols(); ++i)
155  	      (*m_activeColscaleExp)[i] = 0;
156  	
157  	   for(int i = 0; i < lp.nRows(); ++i)
158  	      (*m_activeRowscaleExp)[i] = 0;
159  	
160  	   lp.lp_scaler = this;
161  	}
162  	
163  	
164  	template <class R>
165  	int SPxScaler<R>::computeScaleExp(const SVectorBase<R>& vec,
166  	                                  const DataArray<int>& oldScaleExp) const
167  	{
168  	   R maxi = 0.0;
169  	
170  	   // find largest absolute value after applying existing scaling factors
171  	   for(int i = 0; i < vec.size(); ++i)
172  	   {
173  	      R x = spxAbs(spxLdexp(vec.value(i), oldScaleExp[vec.index(i)]));
174  	
175  	      if(GT(x, maxi, this->tolerances()->epsilon()))
176  	         maxi = x;
177  	   }
178  	
179  	   // empty rows/cols are possible
180  	   if(maxi == 0.0)
181  	      return 0;
182  	   // get exponent corresponding to new scaling factor
183  	   else
184  	   {
185  	      int scaleExp;
186  	      spxFrexp(R(1.0 / maxi), &(scaleExp));
187  	      return scaleExp - 1;
188  	   }
189  	}
190  	
191  	template <class R>
192  	void SPxScaler<R>::applyScaling(SPxLPBase<R>& lp)
193  	{
194  	   assert(lp.nCols() == m_activeColscaleExp->size());
195  	   assert(lp.nRows() == m_activeRowscaleExp->size());
196  	
197  	   DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
198  	   DataArray < int >& rowscaleExp = lp.LPRowSetBase<R>::scaleExp;
199  	
200  	   for(int i = 0; i < lp.nRows(); ++i)
201  	   {
202  	      SVectorBase<R>& vec = lp.rowVector_w(i);
203  	      int exp1;
204  	      int exp2 = rowscaleExp[i];
205  	
206  	      for(int j = 0; j < vec.size(); ++j)
207  	      {
208  	         exp1 = colscaleExp[vec.index(j)];
209  	         vec.value(j) = spxLdexp(vec.value(j), exp1 + exp2);
210  	      }
211  	
212  	      lp.maxRowObj_w(i) = spxLdexp(lp.maxRowObj(i), exp2);
213  	
214  	      if(lp.rhs(i) < R(infinity))
215  	         lp.rhs_w(i) = spxLdexp(lp.rhs_w(i), exp2);
216  	
217  	      if(lp.lhs(i) > R(-infinity))
218  	         lp.lhs_w(i) = spxLdexp(lp.lhs_w(i), exp2);
219  	
220  	      SPxOut::debug(this, "DEBUG: rowscaleExp({}): {}\n", i, exp2);
221  	   }
222  	
223  	   for(int i = 0; i < lp.nCols(); ++i)
224  	   {
225  	      SVectorBase<R>& vec = lp.colVector_w(i);
226  	      int exp1;
227  	      int exp2 = colscaleExp[i];
228  	
229  	      for(int j = 0; j < vec.size(); ++j)
230  	      {
231  	         exp1 = rowscaleExp[vec.index(j)];
232  	         vec.value(j) = spxLdexp(vec.value(j), exp1 + exp2);
233  	      }
234  	
235  	      lp.maxObj_w(i) = spxLdexp(lp.maxObj_w(i), exp2);
236  	
237  	      if(lp.upper(i) < R(infinity))
238  	         lp.upper_w(i) = spxLdexp(lp.upper_w(i), -exp2);
239  	
240  	      if(lp.lower(i) > R(-infinity))
241  	         lp.lower_w(i) = spxLdexp(lp.lower_w(i), -exp2);
242  	
243  	      SPxOut::debug(this, "DEBUG: colscaleExp({}): {}\n", i, exp2);
244  	   }
245  	
246  	   lp.setScalingInfo(true);
247  	   assert(lp.isConsistent());
248  	}
249  	
250  	/// unscale SPxLP
251  	template <class R>
252  	void SPxScaler<R>::unscale(SPxLPBase<R>& lp)
253  	{
254  	   assert(lp.isScaled());
255  	
256  	   const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
257  	   const DataArray < int >& rowscaleExp = lp.LPRowSetBase<R>::scaleExp;
258  	
259  	   for(int i = 0; i < lp.nRows(); ++i)
260  	   {
261  	      SVectorBase<R>& vec = lp.rowVector_w(i);
262  	
263  	      int exp1;
264  	      int exp2 = rowscaleExp[i];
265  	
266  	      for(int j = 0; j < vec.size(); ++j)
267  	      {
268  	         exp1 = colscaleExp[vec.index(j)];
269  	         vec.value(j) = spxLdexp(vec.value(j), -exp1 - exp2);
270  	      }
271  	
272  	      lp.maxRowObj_w(i) = spxLdexp(lp.maxRowObj(i), -exp2);
273  	
274  	      if(lp.rhs(i) < R(infinity))
275  	         lp.rhs_w(i) = spxLdexp(lp.rhs_w(i), -exp2);
276  	
277  	      if(lp.lhs(i) > R(-infinity))
278  	         lp.lhs_w(i) = spxLdexp(lp.lhs_w(i), -exp2);
279  	   }
280  	
281  	   for(int i = 0; i < lp.nCols(); ++i)
282  	   {
283  	      SVectorBase<R>& vec = lp.colVector_w(i);
284  	
285  	      int exp1;
286  	      int exp2 = colscaleExp[i];
287  	
288  	      for(int j = 0; j < vec.size(); ++j)
289  	      {
290  	         exp1 = rowscaleExp[vec.index(j)];
291  	         vec.value(j) = spxLdexp(vec.value(j), -exp1 - exp2);
292  	      }
293  	
294  	      lp.maxObj_w(i) = spxLdexp(lp.maxObj_w(i), -exp2);
295  	
296  	      if(lp.upper(i) < R(infinity))
297  	         lp.upper_w(i) = spxLdexp(lp.upper_w(i), exp2);
298  	
299  	      if(lp.lower(i) > R(-infinity))
300  	         lp.lower_w(i) = spxLdexp(lp.lower_w(i), exp2);
301  	   }
302  	
303  	   lp._isScaled = false;
304  	   assert(lp.isConsistent());
305  	}
306  	
307  	/// returns scaling factor for column \p i
308  	/// todo pass the LP?!
309  	template <class R>
310  	int SPxScaler<R>::getColScaleExp(int i) const
311  	{
312  	   return (*m_activeColscaleExp)[i];
313  	}
314  	
315  	/// returns scaling factor for row \p i
316  	/// todo pass the LP?!
317  	template <class R>
318  	int SPxScaler<R>::getRowScaleExp(int i) const
319  	{
320  	   return (*m_activeRowscaleExp)[i];
321  	}
322  	
323  	/// gets unscaled column \p i
324  	template <class R>
325  	void SPxScaler<R>::getColUnscaled(const SPxLPBase<R>& lp, int i, DSVectorBase<R>& vec) const
326  	{
327  	   assert(lp.isScaled());
328  	   assert(i < lp.nCols());
329  	   assert(i >= 0);
330  	   const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
331  	   const DataArray < int >& rowscaleExp = lp.LPRowSetBase<R>::scaleExp;
332  	
333  	   vec = lp.LPColSetBase<R>::colVector(i);
334  	
335  	   int exp1;
336  	   int exp2 = colscaleExp[i];
337  	
338  	   const SVectorBase<R>& col = lp.colVector(i);
339  	   vec.setMax(col.size());
340  	   vec.clear();
341  	
342  	   for(int j = 0; j < col.size(); j++)
343  	   {
344  	      exp1 = rowscaleExp[col.index(j)];
345  	      vec.add(col.index(j), spxLdexp(col.value(j), -exp1 - exp2));
346  	   }
347  	}
348  	
349  	/// returns maximum absolute value of unscaled column \p i
350  	template <class R>
351  	R SPxScaler<R>::getColMaxAbsUnscaled(const SPxLPBase<R>& lp, int i) const
352  	{
353  	   assert(i < lp.nCols());
354  	   assert(i >= 0);
355  	
356  	   DataArray < int >& colscaleExp = *m_activeColscaleExp;
357  	   DataArray < int >& rowscaleExp = *m_activeRowscaleExp;
358  	   const SVectorBase<R>& colVec = lp.LPColSetBase<R>::colVector(i);
359  	
360  	   R max = 0.0;
361  	   int exp1;
362  	   int exp2 = colscaleExp[i];
363  	
364  	   for(int j = 0; j < colVec.size(); j++)
365  	   {
366  	      exp1 = rowscaleExp[colVec.index(j)];
367  	      R abs = spxAbs(spxLdexp(colVec.value(j), -exp1 - exp2));
368  	
369  	      if(abs > max)
370  	         max = abs;
371  	   }
372  	
373  	   return max;
374  	}
375  	
376  	/// returns minimum absolute value of unscaled column \p i
377  	template <class R>
378  	R SPxScaler<R>::getColMinAbsUnscaled(const SPxLPBase<R>& lp, int i) const
379  	{
380  	   assert(i < lp.nCols());
381  	   assert(i >= 0);
382  	
383  	   DataArray < int >& colscaleExp = *m_activeColscaleExp;
384  	   DataArray < int >& rowscaleExp = *m_activeRowscaleExp;
385  	   const SVectorBase<R>& colVec = lp.LPColSetBase<R>::colVector(i);
386  	
387  	   R min = R(infinity);
388  	   int exp1;
389  	   int exp2 = colscaleExp[i];
390  	
391  	   for(int j = 0; j < colVec.size(); j++)
392  	   {
393  	      exp1 = rowscaleExp[colVec.index(j)];
394  	      R abs = spxAbs(spxLdexp(colVec.value(j), -exp1 - exp2));
395  	
396  	      if(abs < min)
397  	         min = abs;
398  	   }
399  	
400  	   return min;
401  	}
402  	
403  	
404  	/// returns unscaled upper bound \p i
405  	template <class R>
406  	R SPxScaler<R>::upperUnscaled(const SPxLPBase<R>& lp, int i) const
407  	{
408  	   assert(lp.isScaled());
409  	   assert(i < lp.nCols());
410  	   assert(i >= 0);
411  	
412  	   if(lp.LPColSetBase<R>::upper(i) < R(infinity))
413  	   {
414  	      const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
415  	      return spxLdexp(lp.LPColSetBase<R>::upper(i), colscaleExp[i]);
416  	   }
417  	   else
418  	      return lp.LPColSetBase<R>::upper(i);
419  	}
420  	
421  	
422  	/// gets unscaled upper bound vector
423  	template <class R>
424  	void SPxScaler<R>::getUpperUnscaled(const SPxLPBase<R>& lp, VectorBase<R>& vec) const
425  	{
426  	   assert(lp.isScaled());
427  	   assert(lp.LPColSetBase<R>::upper().dim() == vec.dim());
428  	
429  	   const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
430  	
431  	   for(int i = 0; i < lp.LPColSetBase<R>::upper().dim(); i++)
432  	      vec[i] = spxLdexp(lp.LPColSetBase<R>::upper()[i], colscaleExp[i]);
433  	}
434  	
435  	
436  	/// returns unscaled upper bound VectorBase<R> of \p lp
437  	template <class R>
438  	R SPxScaler<R>::lowerUnscaled(const SPxLPBase<R>& lp, int i) const
439  	{
440  	   assert(lp.isScaled());
441  	   assert(i < lp.nCols());
442  	   assert(i >= 0);
443  	
444  	   if(lp.LPColSetBase<R>::lower(i) > R(-infinity))
445  	   {
446  	      const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
447  	      return spxLdexp(lp.LPColSetBase<R>::lower(i), colscaleExp[i]);
448  	   }
449  	   else
450  	      return lp.LPColSetBase<R>::lower(i);
451  	}
452  	
453  	
454  	/// returns unscaled lower bound VectorBase<R> of \p lp
455  	template <class R>
456  	void SPxScaler<R>::getLowerUnscaled(const SPxLPBase<R>& lp, VectorBase<R>& vec) const
457  	{
458  	   assert(lp.isScaled());
459  	   assert(lp.LPColSetBase<R>::lower().dim() == vec.dim());
460  	
461  	   const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
462  	
463  	   for(int i = 0; i < lp.LPColSetBase<R>::lower().dim(); i++)
464  	      vec[i] = spxLdexp(lp.LPColSetBase<R>::lower()[i], colscaleExp[i]);
465  	}
466  	
467  	/// returns unscaled objective function coefficient of \p i
468  	template <class R>
469  	R SPxScaler<R>::maxObjUnscaled(const SPxLPBase<R>& lp, int i) const
470  	{
471  	   assert(lp.isScaled());
472  	   assert(i < lp.nCols());
473  	   assert(i >= 0);
474  	
475  	   const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
476  	
477  	   return spxLdexp(lp.LPColSetBase<R>::maxObj(i), -colscaleExp[i]);
478  	}
479  	
480  	
481  	/// gets unscaled objective function coefficient of \p i
482  	template <class R>
483  	void SPxScaler<R>::getMaxObjUnscaled(const SPxLPBase<R>& lp, VectorBase<R>& vec) const
484  	{
485  	   assert(lp.isScaled());
486  	   assert(lp.LPColSetBase<R>::maxObj().dim() == vec.dim());
487  	
488  	   const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
489  	
490  	   for(int i = 0; i < lp.LPColSetBase<R>::maxObj().dim(); i++)
491  	      vec[i] = spxLdexp(lp.LPColSetBase<R>::maxObj()[i], -colscaleExp[i]);
492  	}
493  	
494  	/// gets unscaled row \p i
495  	template <class R>
496  	void SPxScaler<R>::getRowUnscaled(const SPxLPBase<R>& lp, int i, DSVectorBase<R>& vec) const
497  	{
498  	   assert(lp.isScaled());
499  	   assert(i < lp.nRows());
500  	   assert(i >= 0);
501  	
502  	   const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
503  	   const DataArray < int >& rowscaleExp = lp.LPRowSetBase<R>::scaleExp;
504  	   int exp1;
505  	   int exp2 = rowscaleExp[i];
506  	
507  	   const SVectorBase<R>& row = lp.rowVector(i);
508  	   vec.setMax(row.size());
509  	   vec.clear();
510  	
511  	   for(int j = 0; j < row.size(); j++)
512  	   {
513  	      exp1 = colscaleExp[row.index(j)];
514  	      vec.add(row.index(j), spxLdexp(row.value(j), -exp1 - exp2));
515  	   }
516  	}
517  	
518  	/// returns maximum absolute value of unscaled row \p i
519  	template <class R>
520  	R SPxScaler<R>::getRowMaxAbsUnscaled(const SPxLPBase<R>& lp, int i) const
521  	{
522  	   assert(i < lp.nRows());
523  	   assert(i >= 0);
524  	   DataArray < int >& colscaleExp = *m_activeColscaleExp;
525  	   DataArray < int >& rowscaleExp = *m_activeRowscaleExp;
526  	   const SVectorBase<R>& rowVec = lp.LPRowSetBase<R>::rowVector(i);
527  	
528  	   R max = 0.0;
529  	
530  	   int exp1;
531  	   int exp2 = rowscaleExp[i];
532  	
533  	   for(int j = 0; j < rowVec.size(); j++)
534  	   {
535  	      exp1 = colscaleExp[rowVec.index(j)];
536  	      R abs = spxAbs(spxLdexp(rowVec.value(j), -exp1 - exp2));
537  	
538  	      if(GT(abs, max, this->tolerances()->epsilon()))
539  	         max = abs;
540  	   }
541  	
542  	   return max;
543  	}
544  	
545  	/// returns minimum absolute value of unscaled row \p i
546  	template <class R>
547  	R SPxScaler<R>::getRowMinAbsUnscaled(const SPxLPBase<R>& lp, int i) const
548  	{
549  	   assert(i < lp.nRows());
550  	   assert(i >= 0);
551  	   DataArray < int >& colscaleExp = *m_activeColscaleExp;
552  	   DataArray < int >& rowscaleExp = *m_activeRowscaleExp;
553  	   const SVectorBase<R>& rowVec = lp.LPRowSetBase<R>::rowVector(i);
554  	
555  	   R min = R(infinity);
556  	
557  	   int exp1;
558  	   int exp2 = rowscaleExp[i];
559  	
560  	   for(int j = 0; j < rowVec.size(); j++)
561  	   {
562  	      exp1 = colscaleExp[rowVec.index(j)];
563  	      R abs = spxAbs(spxLdexp(rowVec.value(j), -exp1 - exp2));
564  	
565  	      if(LT(abs, min, this->tolerances()->epsilon()))
566  	         min = abs;
567  	   }
568  	
569  	   return min;
570  	}
571  	
572  	/// returns unscaled right hand side \p i
573  	template <class R>
574  	R SPxScaler<R>::rhsUnscaled(const SPxLPBase<R>& lp, int i) const
575  	{
576  	   assert(lp.isScaled());
577  	   assert(i < lp.nRows());
578  	   assert(i >= 0);
579  	
580  	   if(lp.LPRowSetBase<R>::rhs(i) < R(infinity))
581  	   {
582  	      const DataArray < int >& rowscaleExp = lp.LPRowSetBase<R>::scaleExp;
583  	      return spxLdexp(lp.LPRowSetBase<R>::rhs(i), -rowscaleExp[i]);
584  	   }
585  	   else
586  	      return lp.LPRowSetBase<R>::rhs(i);
587  	}
588  	
589  	
590  	/// gets unscaled right hand side vector
591  	template <class R>
592  	void SPxScaler<R>::getRhsUnscaled(const SPxLPBase<R>& lp, VectorBase<R>& vec) const
593  	{
594  	   assert(lp.isScaled());
595  	   assert(lp.LPRowSetBase<R>::rhs().dim() == vec.dim());
596  	
597  	   for(int i = 0; i < lp.LPRowSetBase<R>::rhs().dim(); i++)
598  	   {
599  	      const DataArray < int >& rowscaleExp = lp.LPRowSetBase<R>::scaleExp;
600  	      vec[i] = spxLdexp(lp.LPRowSetBase<R>::rhs()[i], -rowscaleExp[i]);
601  	   }
602  	}
603  	
604  	
605  	/// returns unscaled left hand side \p i of \p lp
606  	template <class R>
607  	R SPxScaler<R>::lhsUnscaled(const SPxLPBase<R>& lp, int i) const
608  	{
609  	   assert(lp.isScaled());
610  	   assert(i < lp.nRows());
611  	   assert(i >= 0);
612  	
613  	   if(lp.LPRowSetBase<R>::lhs(i) > R(-infinity))
614  	   {
615  	      const DataArray < int >& rowscaleExp = lp.LPRowSetBase<R>::scaleExp;
616  	      return spxLdexp(lp.LPRowSetBase<R>::lhs(i), -rowscaleExp[i]);
617  	   }
618  	   else
619  	      return lp.LPRowSetBase<R>::lhs(i);
620  	}
621  	
622  	/// returns unscaled left hand side VectorBase<R> of \p lp
623  	template <class R>
624  	void SPxScaler<R>::getLhsUnscaled(const SPxLPBase<R>& lp, VectorBase<R>& vec) const
625  	{
626  	   assert(lp.isScaled());
627  	   assert(lp.LPRowSetBase<R>::lhs().dim() == vec.dim());
628  	
629  	   const DataArray < int >& rowscaleExp = lp.LPRowSetBase<R>::scaleExp;
630  	
631  	   for(int i = 0; i < lp.LPRowSetBase<R>::lhs().dim(); i++)
632  	      vec[i] = spxLdexp(lp.LPRowSetBase<R>::lhs()[i], -rowscaleExp[i]);
633  	}
634  	
635  	/// returns unscaled coefficient of \p lp
636  	template <class R>
637  	R SPxScaler<R>::getCoefUnscaled(const SPxLPBase<R>& lp, int row, int col) const
638  	{
639  	   assert(lp.isScaled());
640  	   assert(row < lp.nRows());
641  	   assert(col < lp.nCols());
642  	
643  	   const DataArray < int >& rowscaleExp = lp.LPRowSetBase<R>::scaleExp;
644  	   const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
645  	
646  	   return spxLdexp(lp.colVector(col)[row], - rowscaleExp[row] - colscaleExp[col]);
647  	}
648  	
649  	template <class R>
650  	void SPxScaler<R>::unscalePrimal(const SPxLPBase<R>& lp, VectorBase<R>& x) const
651  	{
652  	   assert(lp.isScaled());
653  	
654  	   const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
655  	
656  	   assert(x.dim() == colscaleExp.size());
657  	
658  	   for(int j = 0; j < x.dim(); ++j)
659  	      x[j] = spxLdexp(x[j], colscaleExp[j]);
660  	}
661  	
662  	template <class R>
663  	void SPxScaler<R>::unscaleSlacks(const SPxLPBase<R>& lp, VectorBase<R>& s) const
664  	{
665  	   assert(lp.isScaled());
666  	
667  	   const DataArray < int >& rowscaleExp = lp.LPRowSetBase<R>::scaleExp;
668  	
669  	   assert(s.dim() == rowscaleExp.size());
670  	
671  	   for(int i = 0; i < s.dim(); ++i)
672  	      s[i] = spxLdexp(s[i], -rowscaleExp[i]);
673  	}
674  	
675  	template <class R>
676  	void SPxScaler<R>::unscaleDual(const SPxLPBase<R>& lp, VectorBase<R>& pi) const
677  	{
678  	   assert(lp.isScaled());
679  	
680  	   const DataArray < int >& rowscaleExp = lp.LPRowSetBase<R>::scaleExp;
681  	
682  	   assert(pi.dim() == rowscaleExp.size());
683  	
684  	   for(int i = 0; i < pi.dim(); ++i)
685  	      pi[i] = spxLdexp(pi[i], rowscaleExp[i]);
686  	}
687  	
688  	template <class R>
689  	void SPxScaler<R>::unscaleRedCost(const SPxLPBase<R>& lp, VectorBase<R>& r) const
690  	{
691  	   assert(lp.isScaled());
692  	
693  	   const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
694  	
695  	   assert(r.dim() == colscaleExp.size());
696  	
697  	   for(int j = 0; j < r.dim(); ++j)
698  	      r[j] = spxLdexp(r[j], -colscaleExp[j]);
699  	}
700  	
701  	template <class R>
702  	void SPxScaler<R>::unscalePrimalray(const SPxLPBase<R>& lp, VectorBase<R>& ray) const
703  	{
704  	   assert(lp.isScaled());
705  	
706  	   const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
707  	
708  	   assert(ray.dim() == colscaleExp.size());
709  	
710  	   for(int j = 0; j < ray.dim(); ++j)
711  	      ray[j] = spxLdexp(ray[j], colscaleExp[j]);
712  	}
713  	
714  	template <class R>
715  	void SPxScaler<R>::unscaleDualray(const SPxLPBase<R>& lp, VectorBase<R>& ray) const
716  	{
717  	   assert(lp.isScaled());
718  	
719  	   const DataArray < int >& rowscaleExp = lp.LPRowSetBase<R>::scaleExp;
720  	
721  	   assert(ray.dim() == rowscaleExp.size());
722  	
723  	   for(int i = 0; i < ray.dim(); ++i)
724  	      ray[i] = spxLdexp(ray[i], rowscaleExp[i]);
725  	}
726  	
727  	template <class R>
728  	void SPxScaler<R>::scaleObj(const SPxLPBase<R>& lp, VectorBase<R>& origObj) const
729  	{
730  	   assert(lp.isScaled());
731  	
732  	   const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
733  	
734  	   for(int i = 0; i < origObj.dim(); ++i)
735  	   {
736  	      origObj[i] = spxLdexp(origObj[i], colscaleExp[i]);
737  	   }
738  	}
739  	
740  	template <class R>
741  	R SPxScaler<R>::scaleObj(const SPxLPBase<R>& lp, int i, R origObj) const
742  	{
743  	   assert(lp.isScaled());
744  	   assert(i < lp.nCols());
745  	   assert(i >= 0);
746  	
747  	   const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
748  	   int exp = colscaleExp[i];
749  	
750  	   return spxLdexp(origObj, exp);
751  	}
752  	
753  	template <class R>
754  	R SPxScaler<R>::scaleElement(const SPxLPBase<R>& lp, int row, int col, R val) const
755  	{
756  	   assert(lp.isScaled());
757  	   assert(col < lp.nCols());
758  	   assert(col >= 0);
759  	   assert(row < lp.nRows());
760  	   assert(row >= 0);
761  	
762  	   const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
763  	   const DataArray < int >& rowscaleExp = lp.LPRowSetBase<R>::scaleExp;
764  	
765  	   return spxLdexp(val, colscaleExp[col] + rowscaleExp[row]);
766  	}
767  	
768  	template <class R>
769  	R SPxScaler<R>::scaleLower(const SPxLPBase<R>& lp, int col, R lower) const
770  	{
771  	   assert(lp.isScaled());
772  	   assert(col < lp.nCols());
773  	   assert(col >= 0);
774  	
775  	   const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
776  	
777  	   return spxLdexp(lower, -colscaleExp[col]);
778  	}
779  	
780  	template <class R>
781  	R SPxScaler<R>::scaleUpper(const SPxLPBase<R>& lp, int col, R upper) const
782  	{
783  	   assert(lp.isScaled());
784  	   assert(col < lp.nCols());
785  	   assert(col >= 0);
786  	
787  	   const DataArray < int >& colscaleExp = lp.LPColSetBase<R>::scaleExp;
788  	
789  	   return spxLdexp(upper, -colscaleExp[col]);
790  	}
791  	
792  	template <class R>
793  	R SPxScaler<R>::scaleLhs(const SPxLPBase<R>& lp, int row, R lhs) const
794  	{
795  	   assert(lp.isScaled());
796  	   assert(row < lp.nRows());
797  	   assert(row >= 0);
798  	
799  	   const DataArray < int >& rowscaleExp = lp.LPRowSetBase<R>::scaleExp;
800  	
801  	   return spxLdexp(lhs, rowscaleExp[row]);
802  	}
803  	
804  	template <class R>
805  	R SPxScaler<R>::scaleRhs(const SPxLPBase<R>& lp, int row, R rhs) const
806  	{
807  	   assert(lp.isScaled());
808  	   assert(row < lp.nRows());
809  	   assert(row >= 0);
810  	
811  	   const DataArray < int >& rowscaleExp = lp.LPRowSetBase<R>::scaleExp;
812  	
813  	   return spxLdexp(rhs, rowscaleExp[row]);
814  	}
815  	
816  	template <class R>
817  	R SPxScaler<R>::minAbsColscale() const
818  	{
819  	   const DataArray < int >& colscaleExp = *m_activeColscaleExp;
820  	
821  	   R mini = R(infinity);
822  	
823  	   for(int i = 0; i < colscaleExp.size(); ++i)
824  	      if(spxAbs(spxLdexp(1.0, colscaleExp[i])) < mini)
825  	         mini = spxAbs(spxLdexp(1.0, colscaleExp[i]));
826  	
827  	   return mini;
828  	}
829  	
830  	template <class R>
831  	R SPxScaler<R>::maxAbsColscale() const
832  	{
833  	   const DataArray < int >& colscaleExp = *m_activeColscaleExp;
834  	
835  	   R maxi = 0.0;
836  	
837  	   for(int i = 0; i < colscaleExp.size(); ++i)
838  	      if(spxAbs(spxLdexp(1.0, colscaleExp[i])) > maxi)
839  	         maxi = spxAbs(spxLdexp(1.0, colscaleExp[i]));
840  	
841  	
842  	   return maxi;
843  	}
844  	
845  	template <class R>
846  	R SPxScaler<R>::minAbsRowscale() const
847  	{
848  	   const DataArray < int >& rowscaleExp = *m_activeRowscaleExp;
849  	
850  	   int mini = std::numeric_limits<int>::max();
851  	
852  	   for(int i = 0; i < rowscaleExp.size(); ++i)
853  	      if(rowscaleExp[i] < mini)
854  	         mini = rowscaleExp[i];
855  	
856  	   return spxLdexp(1.0, mini);
857  	}
858  	
859  	template <class R>
860  	R SPxScaler<R>::maxAbsRowscale() const
861  	{
862  	   const DataArray < int >& rowscaleExp = *m_activeRowscaleExp;
863  	
864  	   int maxi = std::numeric_limits<int>::min();
865  	
866  	   for(int i = 0; i < rowscaleExp.size(); ++i)
867  	      if(rowscaleExp[i] > maxi)
868  	         maxi = rowscaleExp[i];
869  	
870  	   return spxLdexp(1.0, maxi);
871  	}
872  	
873  	/** \f$\max_{j\in\mbox{ cols}}
874  	 *   \left(\frac{\max_{i\in\mbox{ rows}}|a_ij|}
875  	 *              {\min_{i\in\mbox{ rows}}|a_ij|}\right)\f$
876  	 */
877  	template <class R>
878  	R SPxScaler<R>::maxColRatio(const SPxLPBase<R>& lp) const
879  	{
880  	
881  	   R pmax = 0.0;
882  	
883  	   for(int i = 0; i < lp.nCols(); ++i)
884  	   {
885  	      const SVectorBase<R>& vec  = lp.colVector(i);
886  	      R           mini = R(infinity);
887  	      R           maxi = 0.0;
888  	
889  	      for(int j = 0; j < vec.size(); ++j)
890  	      {
891  	         R x = spxAbs(vec.value(j));
892  	
893  	         if(isZero(x, this->tolerances()->epsilon()))
894  	            continue;
895  	
896  	         if(x < mini)
897  	            mini = x;
898  	
899  	         if(x > maxi)
900  	            maxi = x;
901  	      }
902  	
903  	      if(mini == R(infinity))
904  	         continue;
905  	
906  	      R p = maxi / mini;
907  	
908  	      if(p > pmax)
909  	         pmax = p;
910  	   }
911  	
912  	   return pmax;
913  	}
914  	
915  	/** \f$\max_{i\in\mbox{ rows}}
916  	 *   \left(\frac{\max_{j\in\mbox{ cols}}|a_ij|}
917  	 *              {\min_{j\in\mbox{ cols}}|a_ij|}\right)\f$
918  	 */
919  	template <class R>
920  	R SPxScaler<R>::maxRowRatio(const SPxLPBase<R>& lp) const
921  	{
922  	
923  	   R pmax = 0.0;
924  	
925  	   for(int i = 0; i < lp.nRows(); ++i)
926  	   {
927  	      const SVectorBase<R>& vec  = lp.rowVector(i);
928  	      R           mini = R(infinity);
929  	      R           maxi = 0.0;
930  	
931  	      for(int j = 0; j < vec.size(); ++j)
932  	      {
933  	         R x = spxAbs(vec.value(j));
934  	
935  	         if(isZero(x, this->tolerances()->epsilon()))
936  	            continue;
937  	
938  	         if(x < mini)
939  	            mini = x;
940  	
941  	         if(x > maxi)
942  	            maxi = x;
943  	      }
944  	
945  	      if(mini == R(infinity))
946  	         continue;
947  	
948  	      R p = maxi / mini;
949  	
950  	      if(p > pmax)
951  	         pmax = p;
952  	   }
953  	
954  	   return pmax;
955  	}
956  	
957  	template <class R>
958  	void SPxScaler<R>::computeExpVec(const std::vector<R>& vec, DataArray<int>& vecExp)
959  	{
960  	   assert(vec.size() == unsigned(vecExp.size()));
961  	
962  	   for(unsigned i = 0; i < vec.size(); ++i)
963  	   {
964  	      spxFrexp(vec[i], &(vecExp[int(i)]));
965  	      vecExp[int(i)] -= 1;
966  	   }
967  	}
968  	
969  	template <class R>
970  	bool SPxScaler<R>::isConsistent() const
971  	{
972  	#ifdef ENABLE_CONSISTENCY_CHECKS
973  	   return m_activeColscaleExp->isConsistent() && m_activeRowscaleExp->isConsistent();
974  	#else
975  	   return true;
976  	#endif
977  	}
978  	
979  	
980  	} // namespace soplex
981