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 vectorbase.h 26 * @brief Dense vector. 27 */ 28 #ifndef _VECTORBASE_H_ 29 #define _VECTORBASE_H_ 30 31 #include <assert.h> 32 #include <string.h> 33 #include <math.h> 34 #include <iostream> 35 #include "vector" 36 #include "algorithm" 37 38 #include "soplex/spxdefines.h" 39 #include "soplex/stablesum.h" 40 #include "soplex/rational.h" 41 42 namespace soplex 43 { 44 template < class R > class SVectorBase; 45 template < class R > class SSVectorBase; 46 47 /**@brief Dense vector. 48 * @ingroup Algebra 49 * 50 * Class VectorBase provides dense linear algebra vectors. Internally, VectorBase wraps std::vector. 51 * 52 * After construction, the values of a VectorBase can be accessed with the subscript operator[](). Safety is provided by 53 * qchecking of array bound when accessing elements with the subscript operator[]() (only when compiled without \c 54 * -DNDEBUG). 55 * 56 * A VectorBase is distinguished from a simple array of %Reals or %Rationals by providing a set of mathematical 57 * operations. 58 * 59 * The following mathematical operations are provided by class VectorBase (VectorBase \p a, \p b; R \p x): 60 * 61 * <TABLE> 62 * <TR><TD>Operation</TD><TD>Description </TD><TD></TD> </TR> 63 * <TR><TD>\c -= </TD><TD>subtraction </TD><TD>\c a \c -= \c b </TD></TR> 64 * <TR><TD>\c += </TD><TD>addition </TD><TD>\c a \c += \c b </TD></TR> 65 * <TR><TD>\c * </TD><TD>scalar product</TD> 66 * <TD>\c x = \c a \c * \c b </TD></TR> 67 * <TR><TD>\c *= </TD><TD>scaling </TD><TD>\c a \c *= \c x </TD></TR> 68 * <TR><TD>maxAbs() </TD><TD>infinity norm </TD> 69 * <TD>\c a.maxAbs() == \f$\|a\|_{\infty}\f$ </TD></TR> 70 * <TR><TD>minAbs() </TD><TD> </TD> 71 * <TD>\c a.minAbs() == \f$\min |a_i|\f$ </TD></TR> 72 * 73 * <TR><TD>length() </TD><TD>euclidian norm</TD> 74 * <TD>\c a.length() == \f$\sqrt{a^2}\f$ </TD></TR> 75 * <TR><TD>length2()</TD><TD>square norm </TD> 76 * <TD>\c a.length2() == \f$a^2\f$ </TD></TR> 77 * <TR><TD>multAdd(\c x,\c b)</TD><TD>add scaled vector</TD> 78 * <TD> \c a += \c x * \c b </TD></TR> 79 * </TABLE> 80 * 81 * When using any of these operations, the vectors involved must be of the same dimension. Also an SVectorBase \c b is 82 * allowed if it does not contain nonzeros with index greater than the dimension of \c a.q 83 */ 84 template < class R > 85 class VectorBase 86 { 87 88 // VectorBase is a friend of VectorBase of different template type. This is so 89 // that we can do conversions. 90 template <typename S> 91 friend class VectorBase; 92 93 94 protected: 95 96 // ------------------------------------------------------------------------------------------------------------------ 97 /**@name Data */ 98 ///@{ 99 100 /// Values of vector. 101 std::vector<R> val; 102 103 ///@} 104 105 public: 106 107 // ------------------------------------------------------------------------------------------------------------------ 108 /**@name Construction and assignment */ 109 ///@{ 110 111 /// Constructor. 112 /** There is no default constructor since the storage for a VectorBase must be provided externally. Storage must be 113 * passed as a memory block val at construction. It must be large enough to fit at least dimen values. 114 */ 115 116 // Default constructor 117 VectorBase() 118 { 119 // Default constructor 120 ; 121 } 122 123 // Construct from pointer, copies the values into the VectorBase 124 VectorBase(int dimen, R* p_val) 125 { 126 val.assign(p_val, p_val + dimen); 127 } 128 129 // do not convert int to empty vectorbase 130 explicit VectorBase(int p_dimen) 131 { 132 val.resize(p_dimen); 133 } 134 135 // Constructing an element (usually involving casting Real to Rational and 136 // vice versa.) 137 template <typename S> 138 VectorBase(const VectorBase<S>& vec) 139 { 140 this->operator=(vec); 141 } 142 143 // The move constructor 144 VectorBase(const VectorBase<R>&& vec)noexcept: val(std::move(vec.val)) 145 { 146 } 147 148 // Copy constructor 149 VectorBase(const VectorBase<R>& vec): val(vec.val) 150 { 151 } 152 153 154 /// Assignment operator. 155 // Supports assignment from a Rational vector to Real and vice versa 156 template < class S > 157 VectorBase<R>& operator=(const VectorBase<S>& vec) 158 { 159 if((VectorBase<S>*)this != &vec) 160 { 161 val.clear(); 162 val.reserve(vec.dim()); 163 164 for(auto& v : vec.val) 165 { 166 val.push_back(R(v)); 167 } 168 } 169 170 return *this; 171 } 172 173 /// Assignment operator. 174 VectorBase<R>& operator=(const VectorBase<R>& vec) 175 { 176 if(this != &vec) 177 { 178 val.reserve(vec.dim()); 179 180 val = vec.val; 181 } 182 183 return *this; 184 } 185 186 /// Move assignment operator 187 VectorBase<R>& operator=(const VectorBase<R>&& vec) 188 { 189 val = std::move(vec.val); 190 return *this; 191 } 192 193 /// scale and assign 194 VectorBase<R>& scaleAssign(int scaleExp, const VectorBase<R>& vec) 195 { 196 if(this != &vec) 197 { 198 assert(dim() == vec.dim()); 199 200 auto dimen = dim(); 201 202 for(decltype(dimen) i = 0 ; i < dimen; i++) 203 val[i] = spxLdexp(vec[i], scaleExp); 204 205 } 206 207 return *this; 208 } 209 210 /// scale and assign 211 VectorBase<R>& scaleAssign(const int* scaleExp, const VectorBase<R>& vec, bool negateExp = false) 212 { 213 if(this != &vec) 214 { 215 assert(dim() == vec.dim()); 216 217 if(negateExp) 218 { 219 auto dimen = dim(); 220 221 for(decltype(dimen) i = 0; i < dimen; i++) 222 val[i] = spxLdexp(vec[i], -scaleExp[i]); 223 } 224 else 225 { 226 auto dimen = dim(); 227 228 for(decltype(dimen) i = 0; i < dimen; i++) 229 val[i] = spxLdexp(vec[i], scaleExp[i]); 230 } 231 232 } 233 234 return *this; 235 } 236 237 238 /// Assignment operator. 239 /** Assigning an SVectorBase to a VectorBase using operator=() will set all values to 0 except the nonzeros of \p vec. 240 * This is diffent in method assign(). 241 */ 242 template < class S > 243 VectorBase<R>& operator=(const SVectorBase<S>& vec); 244 245 /// Assignment operator. 246 /** Assigning an SSVectorBase to a VectorBase using operator=() will set all values to 0 except the nonzeros of \p 247 * vec. This is diffent in method assign(). 248 */ 249 /**@todo do we need this also in non-template version, because SSVectorBase can be automatically cast to VectorBase? */ 250 template < class S > 251 VectorBase<R>& operator=(const SSVectorBase<S>& vec); 252 253 /// Assign values of \p vec. 254 /** Assigns all nonzeros of \p vec to the vector. All other values remain unchanged. */ 255 template < class S > 256 VectorBase<R>& assign(const SVectorBase<S>& vec); 257 258 /// Assign values of \p vec. 259 /** Assigns all nonzeros of \p vec to the vector. All other values remain unchanged. */ 260 template < class S > 261 VectorBase<R>& assign(const SSVectorBase<S>& vec); 262 263 ///@} 264 265 // ------------------------------------------------------------------------------------------------------------------ 266 /**@name Access */ 267 ///@{ 268 269 /// Dimension of vector. 270 int dim() const 271 { 272 return int(val.size()); 273 } 274 275 /// Return \p n 'th value by reference. 276 R& operator[](int n) 277 { 278 assert(n >= 0 && n < dim()); 279 return val[n]; 280 } 281 282 /// Return \p n 'th value. 283 const R& operator[](int n) const 284 { 285 assert(n >= 0 && n < dim()); 286 return val[n]; 287 } 288 289 /// Equality operator. 290 friend bool operator==(const VectorBase<R>& vec1, const VectorBase<R>& vec2) 291 { 292 return (vec1.val == vec2.val); 293 } 294 295 /// Return underlying std::vector. 296 const std::vector<R>& vec() 297 { 298 return val; 299 } 300 301 ///@} 302 303 // ------------------------------------------------------------------------------------------------------------------ 304 /**@name Arithmetic operations */ 305 ///@{ 306 307 /// Set vector to contain all-zeros (keeping the same length) 308 void clear() 309 { 310 for(auto& v : val) 311 v = 0; 312 } 313 314 /// Addition. 315 template < class S > 316 VectorBase<R>& operator+=(const VectorBase<S>& vec) 317 { 318 assert(dim() == vec.dim()); 319 320 auto dimen = dim(); 321 322 for(decltype(dimen) i = 0; i < dimen; i++) 323 val[i] += vec[i]; 324 325 return *this; 326 } 327 328 /// Addition. 329 template < class S > 330 VectorBase<R>& operator+=(const SVectorBase<S>& vec); 331 332 /// Addition. 333 template < class S > 334 VectorBase<R>& operator+=(const SSVectorBase<S>& vec); 335 336 /// Subtraction. 337 template < class S > 338 VectorBase<R>& operator-=(const VectorBase<S>& vec) 339 { 340 assert(dim() == vec.dim()); 341 342 auto dimen = dim(); 343 344 for(decltype(dimen) i = 0; i < dimen; i++) 345 val[i] -= vec[i]; 346 347 return *this; 348 } 349 350 /// Subtraction. 351 template < class S > 352 VectorBase<R>& operator-=(const SVectorBase<S>& vec); 353 354 /// Subtraction. 355 template < class S > 356 VectorBase<R>& operator-=(const SSVectorBase<S>& vec); 357 358 /// Scaling. 359 template < class S > 360 VectorBase<R>& operator*=(const S& x) 361 { 362 363 auto dimen = dim(); 364 365 for(decltype(dimen) i = 0; i < dimen; i++) 366 val[i] *= x; 367 368 return *this; 369 } 370 371 /// Division. 372 template < class S > 373 VectorBase<R>& operator/=(const S& x) 374 { 375 assert(x != 0); 376 377 auto dimen = dim(); 378 379 for(decltype(dimen) i = 0; i < dimen; i++) 380 val[i] /= x; 381 382 return *this; 383 } 384 385 /// Inner product. 386 R operator*(const VectorBase<R>& vec) const 387 { 388 StableSum<R> x; 389 390 auto dimen = dim(); 391 392 for(decltype(dimen) i = 0; i < dimen; i++) 393 x += val[i] * vec.val[i]; 394 395 return x; 396 } 397 398 /// Inner product. 399 R operator*(const SVectorBase<R>& vec) const; 400 401 /// Inner product. 402 R operator*(const SSVectorBase<R>& vec) const; 403 404 /// Maximum absolute value, i.e., infinity norm. 405 R maxAbs() const 406 { 407 assert(dim() > 0); 408 409 // A helper function for the std::max_element. Because we compare the absolute value. 410 auto absCmpr = [](R a, R b) 411 { 412 return (spxAbs(a) < spxAbs(b)); 413 }; 414 415 auto maxReference = std::max_element(val.begin(), val.end(), absCmpr); 416 417 R maxi = spxAbs(*maxReference); 418 419 assert(maxi >= 0.0); 420 421 return maxi; 422 } 423 424 /// Minimum absolute value. 425 R minAbs() const 426 { 427 assert(dim() > 0); 428 429 // A helper function for the SOPLEX_MIN_element. Because we compare the absolute value. 430 auto absCmpr = [](R a, R b) 431 { 432 return (spxAbs(a) < spxAbs(b)); 433 }; 434 435 auto minReference = SOPLEX_MIN_element(val.begin(), val.end(), absCmpr); 436 437 R mini = spxAbs(*minReference); 438 439 assert(mini >= 0.0); 440 441 return mini; 442 } 443 444 /// Floating point approximation of euclidian norm (without any approximation guarantee). 445 R length() const 446 { 447 return spxSqrt(length2()); 448 } 449 450 /// Squared norm. 451 R length2() const 452 { 453 return (*this) * (*this); 454 } 455 456 /// Addition of scaled vector. 457 template < class S, class T > 458 VectorBase<R>& multAdd(const S& x, const VectorBase<T>& vec) 459 { 460 assert(vec.dim() == dim()); 461 462 auto dimen = dim(); 463 464 for(decltype(dimen) i = 0; i < dimen; i++) 465 val[i] += x * vec.val[i]; 466 467 return *this; 468 } 469 470 /// Addition of scaled vector. 471 template < class S, class T > 472 VectorBase<R>& multAdd(const S& x, const SVectorBase<T>& vec); 473 474 /// Subtraction of scaled vector. 475 template < class S, class T > 476 VectorBase<R>& multSub(const S& x, const SVectorBase<T>& vec); 477 478 /// Addition of scaled vector. 479 template < class S, class T > 480 VectorBase<R>& multAdd(const S& x, const SSVectorBase<T>& vec); 481 482 ///@} 483 484 // ------------------------------------------------------------------------------------------------------------------ 485 /**@name Utilities */ 486 ///@{ 487 488 /// Conversion to C-style pointer. 489 /** This function serves for using a VectorBase in an C-style function. It returns a pointer to the first value of 490 * the array. 491 * 492 * @todo check whether this non-const c-style access should indeed be public 493 */ 494 R* get_ptr() 495 { 496 return val.data(); 497 } 498 499 /// Conversion to C-style pointer. 500 /** This function serves for using a VectorBase in an C-style function. It returns a pointer to the first value of 501 * the array. 502 */ 503 const R* get_const_ptr() const 504 { 505 return val.data(); 506 } 507 508 // Provides access to the iterators of std::vector<R> val 509 typename std::vector<R>::const_iterator begin() const 510 { 511 return val.begin(); 512 } 513 514 typename std::vector<R>::iterator begin() 515 { 516 return val.begin(); 517 } 518 519 // Provides access to the iterators of std::vector<R> val 520 typename std::vector<R>::const_iterator end() const 521 { 522 return val.end(); 523 } 524 525 typename std::vector<R>::iterator end() 526 { 527 return val.end(); 528 } 529 530 // Functions from VectorBase 531 532 // This used to be VectorBase's way of having std::vector's capacity. This 533 // represents the maximum number of elements the std::vector can have without, 534 // needing any more resizing. Bigger than size, mostly. 535 int memSize() const 536 { 537 return int(val.capacity()); 538 } 539 540 /// Resets \ref soplex::VectorBase "VectorBase"'s dimension to \p newdim. 541 void reDim(int newdim, const bool setZero = true) 542 { 543 if(setZero && newdim > dim()) 544 { 545 // Inserts 0 to the rest of the vectors. 546 // 547 // TODO: Is this important after the change of raw pointers to 548 // std::vector. This is just a waste of operations, I think. 549 val.insert(val.end(), newdim - VectorBase<R>::dim(), 0); 550 } 551 else 552 { 553 val.resize(newdim); 554 } 555 556 } 557 558 559 /// Resets \ref soplex::VectorBase "VectorBase"'s memory size to \p newsize. 560 void reSize(int newsize) 561 { 562 assert(newsize > VectorBase<R>::dim()); 563 564 // Problem: This is not a conventional resize for std::vector. This only 565 // updates the capacity, i.e., by pushing elements to the vector after this, 566 // there will not be any (internal) resizes. 567 val.reserve(newsize); 568 } 569 570 // For operations such as vec1 - vec2 571 const VectorBase<R> operator-(const VectorBase<R>& vec) const 572 { 573 assert(vec.dim() == dim()); 574 VectorBase<R> res; 575 res.val.reserve(dim()); 576 577 auto dimen = dim(); 578 579 for(decltype(dimen) i = 0; i < dimen; i++) 580 { 581 res.val.push_back(val[i] - vec[i]); 582 } 583 584 return res; 585 } 586 587 // Addition 588 const VectorBase<R> operator+(const VectorBase<R>& v) const 589 { 590 assert(v.dim() == dim()); 591 VectorBase<R> res; 592 res.val.reserve(dim()); 593 594 auto dimen = dim(); 595 596 for(decltype(dimen) i = 0; i < dimen; i++) 597 { 598 res.val.push_back(val[i] + v[i]); 599 } 600 601 return res; 602 } 603 604 // The negation operator. e.g. -vec1; 605 friend VectorBase<R> operator-(const VectorBase<R>& vec) 606 { 607 VectorBase<R> res; 608 609 res.val.reserve(vec.dim()); 610 611 for(auto& v : vec.val) 612 { 613 res.val.push_back(-(v)); 614 } 615 616 return res; 617 } 618 619 620 ///@} 621 /// Consistency check. 622 bool isConsistent() const 623 { 624 return true; 625 } 626 627 }; 628 629 /// Inner product. 630 template<> 631 inline 632 Rational VectorBase<Rational>::operator*(const VectorBase<Rational>& vec) const 633 { 634 assert(vec.dim() == dim()); 635 636 if(dim() <= 0 || vec.dim() <= 0) 637 return Rational(); 638 639 Rational x = val[0]; 640 x *= vec.val[0]; 641 642 auto dimen = dim(); 643 644 for(decltype(dimen) i = 1; i < dimen; i++) 645 x += val[i] * vec.val[i]; 646 647 return x; 648 } 649 650 } // namespace soplex 651 #endif // _VECTORBASE_H_ 652