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>&nbsp;</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<Rational>" 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
(5) Event move_ctor: User-defined move constructor.
Also see events: [rule_of_five_violation][copy_ctor][copy_assign][remediation][move_assign]
144  	   VectorBase(const VectorBase<R>&& vec)noexcept: val(std::move(vec.val))
145  	   {
146  	   }
147  	
148  	   // Copy constructor
(2) Event copy_ctor: User-defined copy constructor.
Also see events: [rule_of_five_violation][copy_assign][remediation][move_ctor][move_assign]
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.
(3) Event copy_assign: User-defined copy assignment operator.
Also see events: [rule_of_five_violation][copy_ctor][remediation][move_ctor][move_assign]
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
(6) Event move_assign: User-defined move assignment operator.
Also see events: [rule_of_five_violation][copy_ctor][copy_assign][remediation][move_ctor]
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