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