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 >
(1) Event rule_of_five_violation: |
Class "soplex::VectorBase<double>" has a user definition for at least one special function (copy constructor, move constructor, copy assignment, move assignment, destructor) but not all. If one of these functions requires a user definition then the others likely do as well. |
(4) Event remediation: |
Add user-definition for a destructor. |
Also see events: |
[copy_ctor][copy_assign][move_ctor][move_assign] |
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