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   	#include <iostream>
26   	#include <assert.h>
27   	
28   	#include "soplex.h"
29   	#include "soplex/statistics.h"
30   	
31   	namespace soplex
32   	{
33   	/// check scaling of LP
34   	template <class R>
35   	void SoPlexBase<R>::_checkScaling(SPxLPBase<R>* origLP) const
36   	{
37   	   SPX_MSG_INFO1(spxout, spxout << "DEBUG: checking correctness of scaled LP" << std::endl;)
38   	   assert(_realLP->nCols() == origLP->nCols());
39   	   assert(_realLP->nRows() == origLP->nRows());
40   	   assert(_realLP->isScaled() && !origLP->isScaled());
41   	   bool correct = true;
42   	
43   	   SPX_MSG_INFO1(spxout, spxout << "DEBUG: checking rows..." << std::endl;)
44   	
45   	   for(int i = 0; i < origLP->nRows(); ++i)
46   	   {
47   	      assert(EQ(origLP->lhs(i), _realLP->lhsUnscaled(i), this->thesolver.epsilon()));
48   	      assert(EQ(origLP->rhs(i), _realLP->rhsUnscaled(i), this->thesolver.epsilon()));
49   	
50   	      DSVectorBase<R> row;
51   	      _realLP->getRowVectorUnscaled(i, row);
52   	
53   	      assert(origLP->rowVector(i).size() == row.size());
54   	
55   	      for(int j = 0; j < row.size(); ++j)
56   	      {
57   	         if(NE(row.value(j), origLP->rowVector(i).value(j)))
58   	         {
59   	            SPX_MSG_INFO1(spxout, spxout << "DEBUG: scaling error in row " << i << ", col " << j
60   	                          << ": orig " << origLP->rowVector(i).value(j)
61   	                          << ", unscaled: " << row.value(j) << std::endl;)
62   	            correct = false;
63   	         }
64   	      }
65   	   }
66   	
67   	   SPX_MSG_INFO1(spxout, spxout << "DEBUG: checking cols..." << std::endl;)
68   	
69   	   for(int i = 0; i < origLP->nCols(); ++i)
70   	   {
71   	      assert(EQ(origLP->lower(i), _realLP->lowerUnscaled(i), this->thesolver.epsilon()));
72   	      assert(EQ(origLP->upper(i), _realLP->upperUnscaled(i), this->thesolver.epsilon()));
73   	      assert(EQ(origLP->obj(i), _realLP->objUnscaled(i), this->thesolver.epsilon()));
74   	
75   	      DSVectorBase<R> col;
76   	      _realLP->getColVectorUnscaled(i, col);
77   	
78   	      assert(origLP->colVector(i).size() == col.size());
79   	
80   	      for(int j = 0; j < col.size(); ++j)
81   	      {
82   	         if(NE(col.value(j), origLP->colVector(i).value(j), _solver.tolerances()->floatingPointFeastol()))
83   	         {
84   	            SPX_MSG_INFO1(spxout, spxout << "DEBUG: scaling error in col " << i << ", row " << j
85   	                          << ": orig " << origLP->colVector(i).value(j)
86   	                          << ", unscaled: " << col.value(j) << std::endl;)
87   	            correct = false;
88   	         }
89   	      }
90   	   }
91   	
92   	   if(!correct)
93   	   {
94   	      SPX_MSG_INFO1(spxout, spxout << "DEBUG: scaling check failed" << std::endl;)
95   	   }
96   	
97   	   assert(correct);
98   	}
99   	
100  	
101  	template <class R>
102  	void SoPlexBase<R>::_checkBasisScaling()
103  	{
104  	   if(_status != SPxSolverBase<R>::OPTIMAL)
105  	   {
106  	      SPX_MSG_INFO1(spxout, spxout << "DEBUG: skipping test on non optimal bases\n");
107  	      return;
108  	   }
109  	
110  	   assert(&_solver == _realLP);
111  	   VectorBase<R>** binvcol = 0;
112  	   VectorBase<R>** binvrow = 0;
113  	   int* inds = 0;
114  	   int basisdim = _solver.nRows(); // do all operations with regard to the column basis
115  	   bool colrep = (_solver.rep() == SPxSolverBase<R>::COLUMN);
116  	   spx_alloc(binvcol, basisdim);
117  	   spx_alloc(binvrow, basisdim);
118  	   spx_alloc(inds, basisdim);
119  	
120  	   if(colrep)
121  	   {
122  	      SPX_MSG_INFO1(spxout, spxout << "DEBUG: computing columns of inverted basis matrix\n";)
123  	
124  	      // collect columns of the basis inverse
125  	      for(int i = 0; i < basisdim; ++i)
126  	      {
127  	         binvcol[i] = new VectorBase<R>(basisdim);
128  	         binvcol[i]->clear();
129  	         assert(getBasisInverseColReal(i, binvcol[i]->get_ptr(), 0, 0, true));
130  	      }
131  	   }
132  	
133  	   SPX_MSG_INFO1(spxout, spxout << "DEBUG: computing rows of inverted basis matrix\n";)
134  	
135  	   // collect rows of the basis inverse
136  	   for(int i = 0; i < basisdim; ++i)
137  	   {
138  	      binvrow[i] = new VectorBase<R>(basisdim);
139  	      binvrow[i]->clear();
140  	      assert(getBasisInverseRowReal(i, binvrow[i]->get_ptr(), 0, 0, true));
141  	   }
142  	
143  	   if(colrep)
144  	   {
145  	      SPX_MSG_INFO1(spxout, spxout <<
146  	                    "DEBUG: checking columns for identity after multiplying with basis matrix\n";)
147  	
148  	      // multiply with (unscaled) basis matrix and check result (should be unitvecs)
149  	      for(int i = 0; i < basisdim; ++i)
150  	      {
151  	         VectorBase<R> result(*binvcol[i]);
152  	         assert(multBasis(result.get_ptr(), true));
153  	         R sumerror = 0.0;
154  	
155  	         for(int j = 0; j < basisdim; ++j)
156  	         {
157  	            R error = 0.0;
158  	
159  	            if(j != i)
160  	               error = spxAbs(result[j]);
161  	            else
162  	               error = spxAbs(result[j] - 1.0);
163  	
164  	            if(error > _solver.tolerances()->floatingPointFeastol())
165  	               SPX_MSG_INFO1(spxout, spxout << "ERROR: col " << i << " " << j << ", " << result[j] << std::endl);
166  	
167  	            sumerror += error;
168  	         }
169  	
170  	         assert(_solver.rep() == SPxSolverBase<R>::ROW
171  	                || sumerror < _solver.tolerances()->floatingPointFeastol());
172  	      }
173  	   }
174  	
175  	   SPX_MSG_INFO1(spxout, spxout <<
176  	                 "DEBUG: checking rows for identity after multiplying with basis matrix\n";)
177  	
178  	   for(int i = 0; i < basisdim; ++i)
179  	   {
180  	      VectorBase<R> result(*binvrow[i]);
181  	      assert(multBasisTranspose(result.get_ptr(), true));
182  	      R sumerror = 0.0;
183  	
184  	      for(int j = 0; j < basisdim; ++j)
185  	      {
186  	         R error = 0.0;
187  	
188  	         if(j != i)
189  	            error = spxAbs(result[j]);
190  	         else
191  	            error = spxAbs(result[j] - 1.0);
192  	
193  	         if(error > _solver.tolerances()->floatingPointFeastol())
194  	            SPX_MSG_INFO1(spxout, spxout << "ERROR: row " << i << " " << j << ", " << result[j] << std::endl);
195  	
196  	         sumerror += error;
197  	      }
198  	
199  	      assert(sumerror < _solver.tolerances()->floatingPointFeastol());
200  	   }
201  	
202  	   if(_solver.isScaled())
203  	   {
204  	      SPX_MSG_INFO1(spxout, spxout << "DEBUG: unscaling LP\n";)
205  	      //         _solver.setRep(SPxSolverBase<R>::COLUMN);
206  	      _solver.unscaleLPandReloadBasis();
207  	      //         _solver.setBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr());
208  	      //         _solver.solve();
209  	
210  	      VectorBase<R>** binvcol2 = 0;
211  	      VectorBase<R>** binvrow2 = 0;
212  	      spx_alloc(binvcol2, basisdim);
213  	      spx_alloc(binvrow2, basisdim);
214  	
215  	      if(colrep)
216  	      {
217  	         SPX_MSG_INFO1(spxout, spxout << "DEBUG: computing columns of inverted basis matrix again\n";)
218  	
219  	         // collect columns of the basis inverse
220  	         for(int i = 0; i < basisdim; ++i)
221  	         {
222  	            binvcol2[i] = new VectorBase<R>(basisdim);
223  	            binvcol2[i]->clear();
224  	            assert(getBasisInverseColReal(i, binvcol2[i]->get_ptr(), 0, 0, false));
225  	         }
226  	      }
227  	
228  	      SPX_MSG_INFO1(spxout, spxout << "DEBUG: computing rows of inverted basis matrix again\n";)
229  	
230  	      // collect rows of the basis inverse
231  	      for(int i = 0; i < basisdim; ++i)
232  	      {
233  	         binvrow2[i] = new VectorBase<R>(basisdim);
234  	         binvrow2[i]->clear();
235  	         assert(getBasisInverseRowReal(i, binvrow2[i]->get_ptr(), 0, 0, false));
236  	      }
237  	
238  	      SPX_MSG_INFO1(spxout, spxout <<
239  	                    "DEBUG: checking rows and columns of scaled/unscaled inverted of basis matrix\n";)
240  	
241  	      for(int i = 0; i < basisdim; ++i)
242  	      {
243  	         R sumerror = 0.0;
244  	
245  	         for(int j = 0; j < basisdim; ++j)
246  	         {
247  	            if(colrep)
248  	            {
249  	               if(NE((*binvcol[i])[j], (*binvcol2[i])[j], _solver.tolerances()->floatingPointFeastol()))
250  	               {
251  	                  SPX_MSG_INFO1(spxout, spxout << "ERROR: col " << i << " " << j << ", " << (*binvcol[i])[j] << " " <<
252  	                                (*binvcol2[i])[j] << std::endl);
253  	                  sumerror += spxAbs((*binvcol[i])[j] - (*binvcol2[i])[j]);
254  	               }
255  	            }
256  	
257  	            if(NE((*binvrow[i])[j], (*binvrow2[i])[j], _solver.tolerances()->floatingPointFeastol()))
258  	            {
259  	               SPX_MSG_INFO1(spxout, spxout << "ERROR: row " << i << " " << j << ", " << (*binvrow[i])[j] /
260  	                             (*binvrow2[i])[j] << std::endl);
261  	               sumerror += spxAbs((*binvrow[i])[j] - (*binvrow2[i])[j]);
262  	            }
263  	         }
264  	
265  	         assert(sumerror < _solver.tolerances()->floatingPointFeastol());
266  	      }
267  	
268  	      spx_free(binvcol2);
269  	      spx_free(binvrow2);
270  	   }
271  	
272  	   spx_free(inds);
273  	   spx_free(binvrow);
274  	   spx_free(binvcol);
275  	}
276  	
277  	} // namespace soplex
278