/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ /* */ /* This file is part of the HiGHS linear optimization suite */ /* */ /* Written and engineered 2008-2022 at the University of Edinburgh */ /* */ /* Available as open-source under the MIT License */ /* */ /* Authors: Julian Hall, Ivet Galabova, Leona Gottwald and Michael */ /* Feldmeier */ /* */ /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ /**@file simplex/HVectorBase.cpp * @brief */ #include "util/HVectorBase.h" #include #include #include "lp_data/HConst.h" #include "stdio.h" //Just for temporary printf #include "util/HighsCDouble.h" template void HVectorBase::setup(HighsInt size_) { /* * Initialise an HVector instance */ size = size_; count = 0; index.resize(size); array.assign(size, Real{0}); cwork.assign(size + 6400, 0); // MAX invert iwork.assign(size * 4, 0); packCount = 0; packIndex.resize(size); packValue.resize(size); // Initialise three values that are initialised in clear(), but // weren't originally initialised in setup(). Probably doesn't // matter, since clear() is usually called before a vector is // (re-)used. packFlag = false; synthetic_tick = 0; next = 0; } template void HVectorBase::clear() { /* * Clear an HVector instance */ // Standard HVector to clear HighsInt dense_clear = count < 0 || count > size * 0.3; if (dense_clear) { // Treat the array as full if there are no indices or too many // indices array.assign(size, Real{0}); } else { // Zero according to the indices of (possible) nonzeros for (HighsInt i = 0; i < count; i++) { array[index[i]] = 0; } } this->clearScalars(); } template void HVectorBase::clearScalars() { /* * Clear scalars in an HVector instance */ // Reset the flag to indicate when to pack this->packFlag = false; // Zero the number of stored indices this->count = 0; // Zero the synthetic clock for operations with this vector this->synthetic_tick = 0; // Initialise the next value this->next = 0; } template void HVectorBase::tight() { /* * Zero values in Vector.array that do not exceed kHighsTiny in * magnitude, maintaining index if it is well defined */ HighsInt totalCount = 0; using std::abs; if (count < 0) { for (HighsInt my_index = 0; my_index < array.size(); my_index++) if (abs(array[my_index]) < kHighsTiny) array[my_index] = 0; } else { for (HighsInt i = 0; i < count; i++) { const HighsInt my_index = index[i]; const Real& value = array[my_index]; if (abs(value) >= kHighsTiny) { index[totalCount++] = my_index; } else { array[my_index] = Real{0}; } } count = totalCount; } } template void HVectorBase::pack() { /* * Packing (if packFlag set): Pack values/indices in Vector.array * into packValue/Index */ if (!packFlag) return; packFlag = false; packCount = 0; for (HighsInt i = 0; i < count; i++) { const HighsInt ipack = index[i]; packIndex[packCount] = ipack; packValue[packCount] = array[ipack]; packCount++; } } template void HVectorBase::reIndex() { /* * Possibly determine the indices from scratch by passing through * the array */ // Don't do it if there are relatively few nonzeros if (count >= 0 && count <= size * 0.1) return; count = 0; for (HighsInt i = 0; i < size; i++) if ((double)array[i]) index[count++] = i; } template template void HVectorBase::copy(const HVectorBase* from) { /* * Copy from another HVector structure to this instance * The real type of "from" does not need to be the same, but must be * convertible to this HVector's real type. */ clear(); synthetic_tick = from->synthetic_tick; const HighsInt fromCount = count = from->count; const HighsInt* fromIndex = &from->index[0]; const FromReal* fromArray = &from->array[0]; for (HighsInt i = 0; i < fromCount; i++) { const HighsInt iFrom = fromIndex[i]; const FromReal xFrom = fromArray[iFrom]; index[i] = iFrom; array[iFrom] = Real(xFrom); } } template Real HVectorBase::norm2() const { /* * Compute the squared 2-norm of the vector */ const HighsInt workCount = count; const HighsInt* workIndex = &index[0]; const Real* workArray = &array[0]; Real result = Real{0}; for (HighsInt i = 0; i < workCount; i++) { Real value = workArray[workIndex[i]]; result += value * value; } return result; } template template void HVectorBase::saxpy(const RealPivX pivotX, const HVectorBase* pivot) { /* * Add a multiple pivotX of *pivot into this vector, maintaining * indices of nonzeros but not tracking cancellation. * The real types may all be different but must mix in operations and be * convertible to this HVector's real type. */ HighsInt workCount = count; HighsInt* workIndex = &index[0]; Real* workArray = &array[0]; const HighsInt pivotCount = pivot->count; const HighsInt* pivotIndex = &pivot->index[0]; const RealPiv* pivotArray = &pivot->array[0]; using std::abs; for (HighsInt k = 0; k < pivotCount; k++) { const HighsInt iRow = pivotIndex[k]; const Real x0 = workArray[iRow]; const Real x1 = Real(x0 + pivotX * pivotArray[iRow]); if (x0 == Real{0}) workIndex[workCount++] = iRow; workArray[iRow] = (abs(x1) < kHighsTiny) ? kHighsZero : x1; } count = workCount; } template bool HVectorBase::isEqual(const HVectorBase& v0) { if (this->size != v0.size) return false; if (this->count != v0.count) return false; if (this->index != v0.index) return false; if (this->array != v0.array) return false; // if (this->index.size() != v0.index.size()) return false; // for (HighsInt el = 0; el < (HighsInt)this->index.size(); el++) // if (this->index[el] != v0.index[el]) return false; if (this->synthetic_tick != v0.synthetic_tick) return false; return true; } // explicitly instantiate HVectorBase T=double template class HVectorBase; // explicitly instantiate template member function "copy" for // HVectorBase with the type of the copied HVectorBase being either // double or HighsCDouble template void HVectorBase::copy(const HVectorBase*); template void HVectorBase::copy(const HVectorBase*); // explicitly instantiate template member function "saxpy" for // HVectorBase with all four combinations of types for pivot and pivotX: // (double double), (double HighsCDouble), (HighsCDouble double), (HighsCDouble // HighsCDouble) template void HVectorBase::saxpy(const double, const HVectorBase*); template void HVectorBase::saxpy(const double, const HVectorBase*); template void HVectorBase::saxpy(const HighsCDouble, const HVectorBase*); template void HVectorBase::saxpy(const HighsCDouble, const HVectorBase*); // explicitly instantiate HVectorBase T=HighsCDouble template class HVectorBase; // explicitly instantiate template member function "copy" for // HVectorBase with the type of the copied HVectorBase being // either double or HighsCDouble template void HVectorBase::copy(const HVectorBase*); template void HVectorBase::copy(const HVectorBase*); // explicitly instantiate template member function "saxpy" for // HVectorBase with all four combinations of types for pivot and // pivotX: (double double), (double HighsCDouble), (HighsCDouble double), // (HighsCDouble HighsCDouble) template void HVectorBase::saxpy(const double, const HVectorBase*); template void HVectorBase::saxpy( const double, const HVectorBase*); template void HVectorBase::saxpy(const HighsCDouble, const HVectorBase*); template void HVectorBase::saxpy( const HighsCDouble, const HVectorBase*); #if 0 // Todo: Additionally we could add the two specializations to allow pivotX as // int so that integer constants can be used, but I think we can avoid the // additional code bloat and I changed the single place where this is used to // call with 1.0 instead of 1 template void HVectorBase::saxpy(const int, const HVectorBase*); template void HVectorBase::saxpy( const int, const HVectorBase*); template void HVectorBase::saxpy(const int, const HVectorBase*); template void HVectorBase::saxpy(const int, const HVectorBase*); #endif