inmost_sparse.h 43.7 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1 2 3 4 5 6 7 8
#ifndef INMOST_SPARSE_INCLUDED
#define INMOST_SPARSE_INCLUDED


#include "inmost_common.h"

namespace INMOST
{
9 10
	namespace Sparse
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
11
#if defined(USE_SOLVER) || defined(USE_AUTODIFF)
12
		/// Retrive MPI type for row entry type.
Kirill Terekhov's avatar
Kirill Terekhov committed
13
		INMOST_MPI_Type GetRowEntryType();
14
		/// Create MPI type for row entry type.
Kirill Terekhov's avatar
Kirill Terekhov committed
15
		void CreateRowEntryType();
16
		/// Release MPI type.
Kirill Terekhov's avatar
Kirill Terekhov committed
17
		void DestroyRowEntryType();
18
		/// Check whether MPI type was created.
Kirill Terekhov's avatar
Kirill Terekhov committed
19
		bool HaveRowEntryType();
20 21 22 23 24
		/// This class can be used to annotate the matrix.
		/// You can optionally provide a pointer to an object of
		/// this class to Sparse::Matrix::Save or Sparse::HessianMatrix::Save.
		/// Then comments would be added to mtx file next to the
		/// non-empty strings.
Kirill Terekhov's avatar
Kirill Terekhov committed
25 26
		class AnnotationService
		{
27
			interval<INMOST_DATA_ENUM_TYPE,std::string> text; ///< Array of strings corresponding to matrix rows.
Kirill Terekhov's avatar
Kirill Terekhov committed
28
		public:
29 30 31
			/// Create a new annotation for local matrix.
			/// @param start First row of the matrix.
			/// @param end Last row of the matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
32
			AnnotationService(INMOST_DATA_ENUM_TYPE start = 0, INMOST_DATA_ENUM_TYPE end = 0) { if( end != start ) SetInterval(start,end); }
33 34
			/// Copy annotations for the matrix.
			/// @param other Another class for annotation.
Kirill Terekhov's avatar
Kirill Terekhov committed
35
			AnnotationService(const AnnotationService & other) : text(other.text) { }
36 37
			/// Assign annotations for the matrix.
			/// @param other Another class for annotation.
Kirill Terekhov's avatar
Kirill Terekhov committed
38
			AnnotationService & operator = (AnnotationService const & other) {text = other.text; return *this;}
39
			/// Destroy annotations for the matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
40 41 42 43 44
			~AnnotationService() {}
			/// Get the first row index of the distributed matrix interval.
			INMOST_DATA_ENUM_TYPE  GetFirstIndex() const {return text.get_interval_beg();}
			/// Get the last row index of the distributed matrix interval.
			INMOST_DATA_ENUM_TYPE  GetLastIndex() const {return text.get_interval_end();}
45
			/// Check whether the interval of the local matrix was specified and non-empty.
Kirill Terekhov's avatar
Kirill Terekhov committed
46
			bool                   Empty() const {return text.empty();}
47 48 49
			/// Specify interval of the local matrix.
			/// @param beg The first row of the local matrix.
			/// @param end The last row of the local matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
50
			void                   SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end) { text.set_interval_beg(beg); text.set_interval_end(end); }
51 52 53
			/// Retrive interval of the local matrix.
			/// @param start Record the first row of the local matrix into this variable.
			/// @param end Record the last row of the local matrix into this variable.
Kirill Terekhov's avatar
Kirill Terekhov committed
54
			void                   GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const {start = text.get_interval_beg(); end = text.get_interval_end();}
55 56
			/// Retrive the text corresponding to a certain row of the matrix.
			/// @param row Row of the matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
57
			std::string &          GetAnnotation(INMOST_DATA_ENUM_TYPE row) {assert(!text.empty()); return text[row];}
58 59
			/// Retrive the text corresponding to a certain row of the matrix without right of modification.
			/// @param row Row of the matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
60
			const std::string &    GetAnnotation(INMOST_DATA_ENUM_TYPE row) const {assert(!text.empty()); return text[row];}
61 62 63
			/// Specify the text to a certain row of the matrix.
			/// @param row Row of the matrix.
			/// @param str Annotation for the row of the matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
64 65 66 67
			void                   SetAnnotation(INMOST_DATA_ENUM_TYPE row, std::string str) {assert(!text.empty()); text[row] = str;}
		};
#endif
#if defined(USE_SOLVER)
68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101
		/// Distributed vector class.
		/// This class can be used to store both local and distributed dense data of real type.
		/// For example, to form the right-hand side or initial guess to the solution.
		/// @see Solve
		class Vector
		{
		public:
			typedef interval<INMOST_DATA_ENUM_TYPE,INMOST_DATA_REAL_TYPE> Entries; ///< Type for representation of the local array of values.
			typedef Entries::iterator iterator; ///< Type for the iterator over vector of the values.
			typedef Entries::const_iterator const_iterator; ///< Type for the iterator over vector of the constant values.
		private:
			INMOST_MPI_Comm comm; ///< Communicator corresponding to the Vector. Do we need this?
			Entries data; ///< Vector of the vector values on the current processor.
			std::string name; ///< Name of the vector, may be used to specify some parameters.
			bool is_parallel; ///< A flag indicating that the vector has extended range of values via OrderInfo class.
		public:
			/// Main constructor of the Vector class.
			/// @param _name Name of the vector, empty string by default.
			/// @param start Start of the local data interval.
			/// @param end End of the local data interval.
			/// @param _comm Communicator for parallel data exchanges, MPI_COMM_WORLD by default.
			Vector(std::string _name = "", INMOST_DATA_ENUM_TYPE start = 0, INMOST_DATA_ENUM_TYPE end = 0, INMOST_MPI_Comm _comm = INMOST_MPI_COMM_WORLD);
			/// Copy constructor.
			/// @param other Another vector.
			Vector(const Vector & other);
			/// Assignment operator.
			/// @param other Another vector.
			Vector & operator =(Vector const & other);
			/// Delete data of the vector.
			~Vector();
			/// Return reference to i-th element of the vector.
			INMOST_DATA_REAL_TYPE & operator [](INMOST_DATA_ENUM_TYPE i) {return data[i];}
			/// Return i-th element of the vector.
			INMOST_DATA_REAL_TYPE   operator [](INMOST_DATA_ENUM_TYPE i) const {return data[i];}
102 103
			/// Return a block of elements.
			INMOST::Matrix<INMOST_DATA_REAL_TYPE> operator [](const INMOST::AbstractMatrix<INMOST_DATA_INTEGER_TYPE> & rows) const;
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148
			/// Return the global size of the vector.
			INMOST_DATA_ENUM_TYPE  Size() const { return static_cast<INMOST_DATA_ENUM_TYPE>(data.size()); }
			/// Iterator pointing to the first value of the vector.
			iterator             Begin() {return data.begin();}
			/// Iterator pointing to the first constant value of the vector.
			const_iterator       Begin() const {return data.begin();}
			/// Iterator pointing behind the last value of the vector.
			iterator             End() {return data.end();}
			/// Iterator pointing behind the last constant value of the vector.
			const_iterator       End() const {return data.end();}
			/// Test is there any data in the vector.
			bool                 Empty() const {return data.empty();}
			/// Set the start and the end of the distributed vector interval.
			/// @param start The first index of the local part of the vector.
			/// @param end The last index of the local part of the vector.
			void                 SetInterval(INMOST_DATA_ENUM_TYPE   start, INMOST_DATA_ENUM_TYPE   end)       {assert(start<=end); data.set_interval_beg(start); data.set_interval_end(end);}
			/// Get the start and the end of the distributed vector interval.
			/// @param start Record the first index of the local part of the vector into this variable.
			/// @param end Record the last index of the local part of the vector into this variable.
			void                 GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const {start = data.get_interval_beg(); end = data.get_interval_end();}
			/// Move starting position of local indexes
			/// @param shift Number of positions to shift indices.
			void                 ShiftInterval(INMOST_DATA_ENUM_TYPE shift) {data.shift_interval(shift);}
			/// Get the first index of the distributed vector interval.
			INMOST_DATA_ENUM_TYPE  GetFirstIndex() const {return data.get_interval_beg();}
			/// Get the last index of the distributed vector interval.
			INMOST_DATA_ENUM_TYPE  GetLastIndex() const {return data.get_interval_end();}
			/// Get the communicator which the vector is associated with.
			INMOST_MPI_Comm        GetCommunicator() const {return comm;}
			/// Exchange all the data with another vector.
			/// @param other Another vector.
			void Swap(Vector & other) {data.swap(other.data); name.swap(other.name); std::swap(is_parallel,other.is_parallel); std::swap(comm,other.comm);}
			/// Save the distributed vector to a single data file using parallel MPI I/O.
			void                 Save(std::string file);
			/// Load the vector from a single data file using the specified interval.
			/// If interval is not specified, then it will be automatically constructed,
			/// with the about equal block size (the last block may has larger dimension).
			void                 Load(std::string file, INMOST_DATA_ENUM_TYPE mbeg = ENUMUNDEF, INMOST_DATA_ENUM_TYPE mend = ENUMUNDEF, std::string file_ord = "");
			/// Test whether the vector was assigned an extended range of values via OrderInfo class.
			bool                 & isParallel() {return is_parallel;}
			/// Get the vector name specified in the main constructor.
			std::string          GetName() {return name;}
			/// Clear all data of the current vector.
			void Clear() {data.clear();}
		};
Kirill Terekhov's avatar
Kirill Terekhov committed
149
		
150
#endif //defined(USE_SOLVER)
Kirill Terekhov's avatar
Kirill Terekhov committed
151
#if defined(USE_SOLVER) || defined(USE_AUTODIFF)
152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277
		/// Class to store the sparse matrix row.
		/// Represents a sparse vector and facilitates access to individual entries.
		/// This class is used in expressions of automatic differentiation.
		class Row
		{
		public:
			/// Entry of the sparse matrix row.
			typedef struct entry_s
			{
				INMOST_DATA_ENUM_TYPE first;  ///< the column number of the row element.
				INMOST_DATA_REAL_TYPE second; ///< the real value of the row element.
				/// Comparison operator that helps sorting entries.
				bool operator < (const entry_s & other) const { return first < other.first || (first == other.first && second < other.second); }
			} entry;
			/// Assemble an entry of entry_s type.
			/// @param ind Index.
			/// @param val Value.
			__INLINE static entry make_entry(INMOST_DATA_ENUM_TYPE ind, INMOST_DATA_REAL_TYPE val)
			{
				entry ret;
				ret.first = ind;
				ret.second = val;
				return ret;
			}
		private:
			//typedef dynarray<entry,16> Entries; //replace later with more memory-efficient chunk_array, with first chunk in stack
			typedef array<entry> Entries; ///< Container type for 
		public:
			typedef Entries::iterator iterator; ///< Iterator over pairs of index and value.
			typedef Entries::const_iterator const_iterator; ///< Iterator over constant pairs of index and value.
			typedef Entries::reverse_iterator reverse_iterator; ///< Iterator over pairs of index and value running in backward direction.
			typedef Entries::const_reverse_iterator const_reverse_iterator; ///< Iterator over constant pairs of index and value running in backward direction.
		private:
			Entries data; ///< Array of paris of index and value.
		public:
			/// Construct an empty row.
			Row() : data() {}
			/// Copy all data from another row.
			/// @param other Another row.
			Row(const Row & other) : data(other.data) {}
			/// Construct a row from array of pairs of indices and values.
			/// @param pbegin Pointer to the first position in array.
			/// @param pend Pointer behind the last position of array.
			Row(entry * pbegin, entry * pend)  :data(pbegin, pend) {}
			/// Release all data.
			~Row() {}
			/// Copy all data from another row.
			/// @param other Another row.
			Row & operator = (Row const & other) { data = other.data; return *this; }
			/// Finds and returns value with specified index. Adds a new entry if index was not found.
			/// \warning
			/// The operator [] should be used to fill the sparse matrix row, but not to routinely access individual elements of the row.
			/// You can use Sparse::Row::GetIndex and Sparse::Row::GetValue for rapid access to individual elements.
			/// @param i Index.
			/// @return Value corresponding to specified index.
			INMOST_DATA_REAL_TYPE & operator [](INMOST_DATA_ENUM_TYPE i)
			{
				for(Entries::size_type it = 0; it < data.size(); ++it)
					if( data[it].first == i ) return data[it].second;
				data.push_back(make_entry(i,0));
				return data.back().second;
			}
			/// Finds and returns value with specified index. Rises exception on debug and returns extra large value on release if
			/// index was not found.
			/// \warning
			/// The operator [] should be used to fill the sparse matrix row, but not to routinely access individual elements of the row.
			/// You can use Sparse::Row::GetIndex and Sparse::Row::GetValue for rapid access to individual elements.
			/// @param i Index.
			/// @return Value corresponding to specified index.
			INMOST_DATA_REAL_TYPE operator[](INMOST_DATA_ENUM_TYPE i) const
			{
				for (Entries::size_type it = 0; it < data.size(); ++it) if (data[it].first == i) return data[it].second;
				//you should not come here
				assert(false);
				return 1.0e20;
			}
			/// Finds and returns value with specified index.
			/// Returns zero if no entry was found.
			/// @param i Index.
			/// @return Value corresponding to specified index.
			INMOST_DATA_REAL_TYPE get_safe(INMOST_DATA_ENUM_TYPE i) const
			{
				for (Entries::size_type it = 0; it < data.size(); ++it) if (data[it].first == i) return data[it].second;
				return 0.0;
			}
			/// Clear all data of the current row.
			void                    Clear() { data.clear(); }
			/// Exchange all the data with another row.
			/// @param other Another row.
			void                    Swap(Row & other) { data.swap(other.data); }
			/// The size of the sparse row, i.e. the total number of nonzero elements.
			INMOST_DATA_ENUM_TYPE   Size() const { return static_cast<INMOST_DATA_ENUM_TYPE>(data.size()); }
			/// Checks are there any nonzero entries in the row.
			bool                    Empty() const { return data.empty(); }
			/// Retrive an index corresponding to certain position in the array of pairs of index and value.
			/// @param k Position in the array of pairs of index and value.
			/// @return Index corresponding to the position in the array.
			INMOST_DATA_ENUM_TYPE & GetIndex(INMOST_DATA_ENUM_TYPE k) {assert(k < data.size()); return (data.begin()+k)->first;}
			/// Retrive a value corresponding to certain position in the array of pairs of index and value.
			/// @param k Position in the array of pairs of index and value.
			/// @return Value corresponding to the position in the array.
			INMOST_DATA_REAL_TYPE & GetValue(INMOST_DATA_ENUM_TYPE k) {assert(k < data.size()); return (data.begin()+k)->second;}
			/// Retrive an index corresponding to certain position in the array of pairs of index and value.
			/// @param k Position in the array of pairs of index and value.
			/// @return Index corresponding to the position in the array.
			INMOST_DATA_ENUM_TYPE   GetIndex(INMOST_DATA_ENUM_TYPE k) const {assert(k < data.size()); return (data.begin()+k)->first;}
			/// Retrive a value corresponding to certain position in the array of pairs of index and value.
			/// @param k Position in the array of pairs of index and value.
			/// @return Value corresponding to the position in the array.
			INMOST_DATA_REAL_TYPE   GetValue(INMOST_DATA_ENUM_TYPE k) const {assert(k < data.size()); return (data.begin()+k)->second;}
			/// An iterator pointing to the first position in the array of pairs of index and value.
			iterator                Begin() {return data.begin();}
			/// An iterator pointing behind the last position in the array of pairs of index and value.
			iterator                End() {return data.end();}
			/// An iterator pointing to the first position in the array of constant pairs of index and value.
			const_iterator          Begin() const {return data.begin();}
			/// An iterator pointing behind the last position in the array of constant pairs of index and value.
			const_iterator          End() const {return data.end();}
			/// An iterator pointing to the last position in the array of pairs of index and value.
			reverse_iterator        rBegin() { return data.rbegin(); }
			/// An iterator pointing before the first position in the array of pairs of index and value.
			reverse_iterator        rEnd() { return data.rend(); }
			/// An iterator pointing to the last position in the array of constant pairs of index and value.
			const_reverse_iterator  rBegin() const { return data.rbegin(); }
			/// An iterator pointing before the first position in the array of constant pairs of index and value.
			const_reverse_iterator  rEnd() const { return data.rend(); }
Kirill Terekhov's avatar
Kirill Terekhov committed
278
#if defined(USE_SOLVER)
279 280
			/// Return the scalar product of the current sparse row by a dense Vector.
			INMOST_DATA_REAL_TYPE   RowVec(Vector & x) const; // returns A(row) * x
Kirill Terekhov's avatar
Kirill Terekhov committed
281
#endif
282 283 284 285 286 287 288 289 290 291 292 293 294 295
			/// An optimized assignment of the row, when the content of the source row may not be preserved.
			/// @param source Source raw where to get the contents.
			void                    MoveRow(Row & source) {data = source.data;} //here move constructor and std::move may be used in future
			/// Set the vector entries by zeroes.
			void                    Zero() {for(iterator it = Begin(); it != End(); ++it) it->second = 0;}
			/// Push specified element into sparse row.
			/// This function should be used only if the index is not repeated in the row.
			void                    Push(INMOST_DATA_ENUM_TYPE ind, INMOST_DATA_REAL_TYPE val) {data.push_back(make_entry(ind,val));}
			/// Resize row to specified size.
			/// It is intended to be used together with non-const Row::GetIndex and Row::GetValue
			/// that allow for the modification of individual entries.
			/// @param size New size of the row.
			void                    Resize(INMOST_DATA_ENUM_TYPE size) {data.resize(size);}
			/// Output all entries of the row.
296
			void                    Print(double eps = -1) const
297
			{
298 299 300
				int k = 0;
				for(const_iterator it = Begin(); it != End(); ++it) if( fabs(it->second) > eps ) {std::cout << "(" << it->first << "," << it->second << ") "; k++; }
				if( k ) std::cout << std::endl;
301 302 303 304 305 306 307 308 309 310 311 312
			}
			/// Check whether the row is sorted.
			bool                    isSorted() const;
			/// Add up two rows. Performs operation output=alpha*left+beta*right.
			/// @param alpha Coefficient to multiply the left row.
			/// @param left The left row.
			/// @param beta Coefficient to multiply the right row.
			/// @param right The right row.
			/// @param output Record result in this vector.
			static void             MergeSortedRows(INMOST_DATA_REAL_TYPE alpha, const Row & left, INMOST_DATA_REAL_TYPE beta, const Row & right, Row & output);
		};
		
Kirill Terekhov's avatar
Kirill Terekhov committed
313 314 315
#endif //defined(USE_SOLVER) || defined(USE_AUTODIFF)
		
#if defined(USE_SOLVER) || defined(USE_AUTODIFF)
316 317 318 319 320 321 322 323 324 325 326 327 328 329 330
		/// Class to store the compressed symmetric matrix of a hessian row.
		class HessianRow
		{
		public:
			/// Entry of the sparse matrix row.
			typedef struct hessian_index_s
			{
				INMOST_DATA_ENUM_TYPE first;
				INMOST_DATA_ENUM_TYPE second;
				bool operator < (const hessian_index_s & other) const {return first < other.first || (first == other.first && second < other.second); }
				bool operator ==(const hessian_index_s & other) const {return first == other.first && second == other.second;}
			} index;
			__INLINE static index make_index(INMOST_DATA_ENUM_TYPE _first, INMOST_DATA_ENUM_TYPE _second)
			{
				index ret;
331 332
				ret.first = _first >_second ? _second : _first;
				ret.second = _first < _second ? _second : _first;
333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427
				return ret;
			}
			typedef struct hessian_entry_s
			{
				index first;  ///< the column number of the row element.
				INMOST_DATA_REAL_TYPE second; ///< the real value of the row element.
				bool operator < (const hessian_entry_s & other) const { return first < other.first || (first == other.first && second < other.second); }
			} entry;
			__INLINE static entry make_entry(index ind, INMOST_DATA_REAL_TYPE val)
			{
				entry ret;
				ret.first = ind;
				ret.second = val;
				return ret;
			}
		private:
			typedef array<entry> Entries; //replace later with more memory-efficient chunk_array, with first chunk in stack
			//typedef dynarray<entry,8> Entries;
		public:
			typedef Entries::iterator iterator;
			typedef Entries::const_iterator const_iterator;
			typedef Entries::reverse_iterator reverse_iterator;
			typedef Entries::const_reverse_iterator const_reverse_iterator;
		private:
			Entries data;
		public:
			HessianRow() : data() {}
			HessianRow(const HessianRow & other) : data(other.data) {}
			HessianRow(entry * pbegin, entry * pend) : data(pbegin,pend) {}
			~HessianRow() {}
			HessianRow & operator = (HessianRow const & other) {data = other.data; return *this;}
			/// The operator [] used to fill the sparse matrix row, but not to access individual elements of the row.
			INMOST_DATA_REAL_TYPE & operator [](index i) // use to fill matrix, not to access individual elements
			{
				for(Entries::size_type it = 0; it < data.size(); ++it)
					if( data[it].first == i ) return data[it].second;
				data.push_back(make_entry(i,0));
				return data.back().second;
			}
			/// The operator [] used to access individual elements of the row.
			INMOST_DATA_REAL_TYPE operator[](index i) const
			{
				for (Entries::size_type it = 0; it < data.size(); ++it) if (data[it].first == i) return data[it].second;
				//you should not come here
				assert(false);
				return 1.0e20;
			}
			INMOST_DATA_REAL_TYPE get_safe(index i) const
			{
				for (Entries::size_type it = 0; it < data.size(); ++it) if (data[it].first == i) return data[it].second;
				return 0.0;
			}
			//void           Reserve(INMOST_DATA_ENUM_TYPE num) { data.reserve(num);}
			/// Clear all data of the current row.
			void                    Clear() { data.clear(); }
			void                    Swap(HessianRow & other) { data.swap(other.data); }
			/// The size of the sparse row, i.e. the total number of nonzero elements.
			INMOST_DATA_ENUM_TYPE   Size() const { return static_cast<INMOST_DATA_ENUM_TYPE>(data.size()); }
			bool                    Empty() const { return data.empty(); }
			index                 & GetIndex(INMOST_DATA_ENUM_TYPE k) {assert(k < data.size()); return (data.begin()+k)->first;}
			INMOST_DATA_REAL_TYPE & GetValue(INMOST_DATA_ENUM_TYPE k) {assert(k < data.size()); return (data.begin()+k)->second;}
			index                   GetIndex(INMOST_DATA_ENUM_TYPE k) const {assert(k < data.size()); return (data.begin()+k)->first;}
			INMOST_DATA_REAL_TYPE   GetValue(INMOST_DATA_ENUM_TYPE k) const {assert(k < data.size()); return (data.begin()+k)->second;}
			iterator                Begin() {return data.begin();}
			iterator                End() {return data.end();}
			const_iterator          Begin() const {return data.begin();}
			const_iterator          End() const {return data.end();}
			reverse_iterator        rBegin() { return data.rbegin(); }
			reverse_iterator        rEnd() { return data.rend(); }
			const_reverse_iterator  rBegin() const { return data.rbegin(); }
			const_reverse_iterator  rEnd() const { return data.rend(); }
			/// Return the scalar product of the current sparse row by a dense Vector.
			void                    RowVec(INMOST_DATA_REAL_TYPE alpha, const Row & rU, INMOST_DATA_REAL_TYPE beta, Row & rJ) const; // returns A(row) * x
			void                    MoveRow(HessianRow & new_pos) {data = new_pos.data;} //here move constructor and std::move may be used in future
			/// Set the vector entries by zeroes.
			void                    Zero() {for(iterator it = Begin(); it != End(); ++it) it->second = 0;}
			/// Push specified element into sparse row.
			/// This function should be used only if the index is not repeated in the row.
			void                    Push(index ind, INMOST_DATA_REAL_TYPE val) {data.push_back(make_entry(ind,val));}
			/// Resize row to specified size.
			/// It is intended to be used together with non-const Row::GetIndex and Row::GetValue
			/// that allow for the modification of individual entries.
			/// @param size New size of the row.
			void                    Resize(INMOST_DATA_ENUM_TYPE size) {data.resize(size);}
			void                    Print() const
			{
				for(const_iterator it = Begin(); it != End(); ++it) std::cout << "(" << it->first.first << "," << it->first.second << "," << it->second << ") ";
				std::cout << std::endl;
			}
			bool                    isSorted() const;
			/// output = alpha * left + beta *right
			static void             MergeSortedRows(INMOST_DATA_REAL_TYPE alpha, const HessianRow & left, INMOST_DATA_REAL_TYPE beta, const HessianRow & right, HessianRow & output);
			static void             MergeJacobianHessian(INMOST_DATA_REAL_TYPE a, const Row & JL, const Row & JR, INMOST_DATA_REAL_TYPE b, const HessianRow & HL, INMOST_DATA_REAL_TYPE c, const HessianRow & HR, HessianRow & output);
			static void             MergeJacobianHessian(INMOST_DATA_REAL_TYPE a, const Row & JL, const Row & JR, INMOST_DATA_REAL_TYPE b, const HessianRow & H, HessianRow & output);
		};
Kirill Terekhov's avatar
Kirill Terekhov committed
428 429 430 431 432
		
#endif //defined(USE_SOLVER) || defined(USE_AUTODIFF)
		
#if defined(USE_SOLVER)
		
433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524
		/// This class can be used for shared access to matrix with OpenMP.
		class LockService
		{
			interval<INMOST_DATA_ENUM_TYPE,INMOST_OMP_LOCK_T> locks;
			void DestroyLocks();
		public:
			LockService(INMOST_DATA_ENUM_TYPE start = 0, INMOST_DATA_ENUM_TYPE end = 0) { if( end != start ) SetInterval(start,end); }
			LockService(const LockService & other) { SetInterval(other.GetFirstIndex(),other.GetLastIndex()); }
			LockService & operator = (LockService const & other) { SetInterval(other.GetFirstIndex(),other.GetLastIndex()); return *this;}
			~LockService() {DestroyLocks();}
			/// Get the first row index of the distributed matrix interval.
			INMOST_DATA_ENUM_TYPE  GetFirstIndex() const;
			/// Get the last row index of the distributed matrix interval.
			INMOST_DATA_ENUM_TYPE  GetLastIndex() const;
			bool Empty() const;
			void SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end);
			void GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const;
			bool HaveLocks() const;
			bool Lock(INMOST_DATA_ENUM_TYPE row);
			bool TestLock(INMOST_DATA_ENUM_TYPE row);
			bool UnLock(INMOST_DATA_ENUM_TYPE row);
		};
		
		/// Class to store the distributed sparse matrix by compressed rows.
		/// The format used to store sparse matrix is analogous to Compressed Row Storage format (CRS).
		/// @see http://netlib.org/linalg/html_templates/node91.html
		class Matrix
		{
			typedef interval<INMOST_DATA_ENUM_TYPE,Row> Rows;
		public:
			typedef Rows::iterator iterator;
			typedef Rows::const_iterator const_iterator;
		private:
			INMOST_MPI_Comm comm;
			Rows data;
			std::string name;
			bool is_parallel;
		public:
			/// Main constructor of the Matrix class.
			/// @param _name Name of the matrix, empty string by default.
			/// @param start Start of the local data interval.
			/// @param end End of the local data interval.
			/// @param _comm Communicator for parallel data exchanges, MPI_COMM_WORLD by default.
			Matrix(std::string _name = "", INMOST_DATA_ENUM_TYPE start = 0, INMOST_DATA_ENUM_TYPE end = 0, INMOST_MPI_Comm _comm = INMOST_MPI_COMM_WORLD);
			Matrix(const Matrix & other);
			Matrix & operator =(Matrix const & other);
			~Matrix();
			/// Return reference to i-th Row of the matrix.
			Row & operator [](INMOST_DATA_ENUM_TYPE i) {return data[i];}
			/// Return reference to i-th Row of the matrix.
			const Row & operator [](INMOST_DATA_ENUM_TYPE i) const {return data[i];}
			/// Return the total number of rows in the matrix.
			INMOST_DATA_ENUM_TYPE  Size() const { return static_cast<INMOST_DATA_ENUM_TYPE>(data.size()); }
			bool                   Empty() const {return data.empty();}
			iterator               Begin() {return data.begin();}
			iterator               End()   {return data.end();}
			const_iterator         Begin() const {return data.begin();}
			const_iterator         End() const   {return data.end();}
			/// Set the start and the end row numbers of the distributed matrix interval.
			void                   SetInterval(INMOST_DATA_ENUM_TYPE   start, INMOST_DATA_ENUM_TYPE   end) {data.set_interval_beg(start); data.set_interval_end(end);}
			/// Get the start and the end row numbers of the distributed matrix interval.
			void                   GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const {start = data.get_interval_beg(); end = data.get_interval_end();}
			void                   ShiftInterval(INMOST_DATA_ENUM_TYPE shift) {data.shift_interval(shift);}
			/// Get the first row index of the distributed matrix interval.
			INMOST_DATA_ENUM_TYPE  GetFirstIndex() const {return data.get_interval_beg();}
			/// Get the last row index of the distributed matrix interval.
			INMOST_DATA_ENUM_TYPE  GetLastIndex() const {return data.get_interval_end();}
			/// Get the communicator which the matrix is associated with.
			INMOST_MPI_Comm        GetCommunicator() const {return comm;}
			void                   MoveRows(INMOST_DATA_ENUM_TYPE from, INMOST_DATA_ENUM_TYPE to, INMOST_DATA_ENUM_TYPE size); //for parallel
			void                   Swap(Matrix & other);
			/// Matrix-vector product of the form: y = alpha*A*x + beta * y.
			/// @param y Input/output vector.
			void MatVec(INMOST_DATA_REAL_TYPE alpha, Vector & x, INMOST_DATA_REAL_TYPE beta, Vector & y) const;
			/// Matrix-vector product with transposed matrix of the form: y = alpha*A^T*x + beta * y.
			/// @param y Input/output vector.
			void MatVecTranspose(INMOST_DATA_REAL_TYPE alpha, Vector & x, INMOST_DATA_REAL_TYPE beta, Vector & y) const;
			/// Clear all data of the matrix.
			void Clear() {for(Matrix::iterator it = Begin(); it != End(); ++it) it->Clear(); data.clear();}
			/// Load the matrix from a single data file in MTX format using the specified interval.
			/// If interval is not specified, then it will be automatically constructed,
			/// with the about equal block size (the last block may has larger dimension).
			void				 Load(std::string file, INMOST_DATA_ENUM_TYPE beg = ENUMUNDEF, INMOST_DATA_ENUM_TYPE end = ENUMUNDEF, std::string file_ord = "");
			/// Save the distributed matrix to a single data file in MTX format using parallel MPI I/O.
			/// @see http://math.nist.gov/MatrixMarket/formats.html
			void                 Save(std::string file, const AnnotationService * annotation = NULL);
			/// Check that matrix is in parallel state.
			bool &               isParallel() { return is_parallel; }
			const bool &         isParallel() const { return is_parallel; }
			/// Get the matrix name specified in the main constructor.
			std::string          GetName() const {return name;}
		};
Kirill Terekhov's avatar
Kirill Terekhov committed
525 526 527 528 529
		
#endif //defined(USE_SOLVER)
		
#if defined(USE_SOLVER)
		
530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596
		/// Class to store the distributed sparse hessian hyper matrix by compressed symmetric matrices.
		class HessianMatrix
		{
		public:
			typedef interval<INMOST_DATA_ENUM_TYPE,HessianRow> HessianRows;
			typedef HessianRows::iterator iterator;
			typedef HessianRows::const_iterator const_iterator;
		private:
			INMOST_MPI_Comm comm;
			HessianRows data;
			std::string name;
			bool is_parallel;
		public:
			/// Main constructor of the Matrix class.
			/// @param _name Name of the matrix, empty string by default.
			/// @param start Start of the local data interval.
			/// @param end End of the local data interval.
			/// @param _comm Communicator for parallel data exchanges, MPI_COMM_WORLD by default.
			HessianMatrix(std::string _name = "", INMOST_DATA_ENUM_TYPE start = 0, INMOST_DATA_ENUM_TYPE end = 0, INMOST_MPI_Comm _comm = INMOST_MPI_COMM_WORLD);
			HessianMatrix(const HessianMatrix & other);
			HessianMatrix & operator =(HessianMatrix const & other);
			~HessianMatrix();
			/// Return reference to i-th Row of the matrix.
			HessianRow & operator [](INMOST_DATA_ENUM_TYPE i) {return data[i];}
			/// Return reference to i-th Row of the matrix.
			const HessianRow & operator [](INMOST_DATA_ENUM_TYPE i) const {return data[i];}
			/// Return the total number of rows in the matrix.
			INMOST_DATA_ENUM_TYPE  Size() const { return static_cast<INMOST_DATA_ENUM_TYPE>(data.size()); }
			bool                   Empty() const {return data.empty();}
			iterator               Begin() {return data.begin();}
			iterator               End()   {return data.end();}
			const_iterator         Begin() const {return data.begin();}
			const_iterator         End() const   {return data.end();}
			/// Set the start and the end row numbers of the distributed matrix interval.
			void                   SetInterval(INMOST_DATA_ENUM_TYPE   start, INMOST_DATA_ENUM_TYPE   end) {data.set_interval_beg(start); data.set_interval_end(end);}
			/// Get the start and the end row numbers of the distributed matrix interval.
			void                   GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const {start = data.get_interval_beg(); end = data.get_interval_end();}
			void                   ShiftInterval(INMOST_DATA_ENUM_TYPE shift) {data.shift_interval(shift);}
			/// Get the first row index of the distributed matrix interval.
			INMOST_DATA_ENUM_TYPE  GetFirstIndex() const {return data.get_interval_beg();}
			/// Get the last row index of the distributed matrix interval.
			INMOST_DATA_ENUM_TYPE  GetLastIndex() const {return data.get_interval_end();}
			/// Get the communicator which the matrix is associated with.
			INMOST_MPI_Comm        GetCommunicator() const {return comm;}
			void                   MoveRows(INMOST_DATA_ENUM_TYPE from, INMOST_DATA_ENUM_TYPE to, INMOST_DATA_ENUM_TYPE size); //for parallel
			void                   Swap(HessianMatrix & other);
			/// HyperMatrix-Matrix product of the form: J = alpha*H*U + beta * J.
			/// J - jacobian sparse matrix
			/// H - hessian symmetric sparse hypermatrix
			/// U - update sparse matrix
			/// @param J Input/output Matrix.
			void MatVec(INMOST_DATA_REAL_TYPE alpha, const Matrix & U, INMOST_DATA_REAL_TYPE beta, Matrix & J) const;
			/// Clear all data of the matrix.
			void Clear() {for(HessianMatrix::iterator it = Begin(); it != End(); ++it) it->Clear(); data.clear();}
			/// Load the matrix from a single data file in MTX format using the specified interval.
			/// If interval is not specified, then it will be automatically constructed,
			/// with the about equal block size (the last block may has larger dimension).
			void				 Load(std::string file, INMOST_DATA_ENUM_TYPE beg = ENUMUNDEF, INMOST_DATA_ENUM_TYPE end = ENUMUNDEF, std::string file_ord = "");
			/// Save the distributed matrix to a single data file in MTX format using parallel MPI I/O.
			/// @see http://math.nist.gov/MatrixMarket/formats.html
			void                 Save(std::string file, const AnnotationService * annotation = NULL);
			/// Check that matrix is in parallel state.
			bool &               isParallel() { return is_parallel; }
			const bool &         isParallel() const { return is_parallel; }
			/// Get the matrix name specified in the main constructor.
			std::string          GetName() const {return name;}
		};
Kirill Terekhov's avatar
Kirill Terekhov committed
597 598 599 600
		
#endif //defined(USE_SOLVER)
		
#if defined(USE_SOLVER) || defined(USE_AUTODIFF)
601
		/// This class may be used to sum multiple sparse rows.
602 603 604 605 606
		/// \warning
		/// In parallel column indices of the matrix may span wider then
		/// local row indices, to prevent any problem you can safely
		/// set total size of the matrix as interval of the RowMerger.
		/// For efficiency you should use RowMerger::SetNonlocal function
Kirill Terekhov's avatar
Kirill Terekhov committed
607
		/// to provide information about non-local elements.
608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639
		class RowMerger
		{
		public:
			static const INMOST_DATA_ENUM_TYPE EOL = ENUMUNDEF-1; ///< End of linked list.
			static const INMOST_DATA_ENUM_TYPE UNDEF = ENUMUNDEF; ///< Value not defined in linked list.
			class iterator
			{
			private:
				INMOST_DATA_ENUM_TYPE pos;
				interval< INMOST_DATA_ENUM_TYPE, Row::entry > * LinkedList; ///< Link to associated storage for linked list.
				iterator(interval< INMOST_DATA_ENUM_TYPE, Row::entry > * pLinkedList) : pos(pLinkedList->get_interval_beg()), LinkedList(pLinkedList) {}
			public:
				iterator(const iterator & other) : pos(other.pos), LinkedList(other.LinkedList) {}
				~iterator() {}
				INMOST_DATA_REAL_TYPE & operator *() {return (*LinkedList)[pos].second;}
				INMOST_DATA_REAL_TYPE operator *() const {return (*LinkedList)[pos].second;}
				INMOST_DATA_REAL_TYPE * operator ->() {return &(*LinkedList)[pos].second;}
				const INMOST_DATA_REAL_TYPE * operator ->() const {return &(*LinkedList)[pos].second;}
				iterator & operator ++(){ pos = (*LinkedList)[pos].first; return *this;}
				iterator operator ++(int){ iterator ret(LinkedList); ret.pos = (*LinkedList)[pos].first; return ret; }
				iterator & operator = (const iterator & other) {LinkedList = other.LinkedList; pos = other.pos; return *this;}
				bool operator ==(const iterator & other) const {return LinkedList == other.LinkedList && pos == other.pos;}
				bool operator !=(const iterator & other) const {return LinkedList != other.LinkedList || pos != other.pos;}
				bool operator < (const iterator & other) const {return LinkedList == other.LinkedList && pos < other.pos;}
				bool operator <=(const iterator & other) const {return LinkedList == other.LinkedList && pos <= other.pos;}
				bool operator > (const iterator & other) const {return LinkedList == other.LinkedList && pos > other.pos;}
				bool operator >=(const iterator & other) const {return LinkedList == other.LinkedList && pos >= other.pos;}
				friend class RowMerger;
			};
		private:
			bool Sorted; ///< Contents of linked list should be sorted.
			INMOST_DATA_ENUM_TYPE Nonzeros; ///< Number of nonzero in linked list.
640 641
			INMOST_DATA_ENUM_TYPE IntervalBeg; ///< Begin of global interval of owned index interval
			INMOST_DATA_ENUM_TYPE IntervalEnd;  ///< End of global interval of owned index interval
642
			interval< INMOST_DATA_ENUM_TYPE, Row::entry > LinkedList; ///< Storage for linked list.
643 644
			std::vector< INMOST_DATA_ENUM_TYPE > NonlocalPre; ///< List of global indices, that are to the left of owned index interval
			std::vector< INMOST_DATA_ENUM_TYPE > NonlocalPost; ///< List of global indices, that are to the right of owned index interval
645
			std::map< INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE > Nonlocal; ///< Additional space
646
		public:
647 648 649 650 651 652 653 654 655 656 657 658 659 660
			/// This function converts global index into local index.
			/// @param pos Global index.
			/// @return Local index.
			INMOST_DATA_ENUM_TYPE MapIndex(INMOST_DATA_ENUM_TYPE pos) const;
			/// This function converts local index into global index.
			/// @param pos Local index.
			/// @return Global index.
			INMOST_DATA_ENUM_TYPE UnmapIndex(INMOST_DATA_ENUM_TYPE pos) const;
			/// This function provides information about additional non-local indices.
			/// \warning
			/// All contents of linked list will be lost.
			/// @param Pre Non-local indices that go before IntervalBegin.
			/// @param Post Non-local indices that follow IntervalEnd.
			void SetNonlocal(const std::vector<INMOST_DATA_ENUM_TYPE> & Pre, const std::vector<INMOST_DATA_ENUM_TYPE> & Post);
661 662 663 664 665 666 667
			/// Default constructor without size specified.
			RowMerger();
			/// Constructor with size specified.
			/// @param interval_begin First index in linked list.
			/// @param interval_end Last index in linked list.
			/// @param Sorted Result should be sorted or not.
			RowMerger(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool Sorted = true);
668 669 670 671 672 673 674 675
			/// Constructor with size and non-local mapping specified.
			/// @param interval_begin First index in linked list.
			/// @param interval_end Last index in linked list.
			/// @param Pre Nonlocal indices before First index in linked list.
			/// @param Post Nonlocal indices after Last index in linked list.
			/// @param Sorted Result should be sorted or not.
			RowMerger(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, const std::vector<INMOST_DATA_ENUM_TYPE> & Pre, const std::vector<INMOST_DATA_ENUM_TYPE> & Post, bool Sorted = true);
			/// Destructor.
676 677 678 679 680 681 682 683
			~RowMerger();
			/// Resize linked list for new interval.
			/// \warning
			/// All contents of linked list will be lost after resize.
			/// @param interval_begin First index in linked list.
			/// @param interval_end Last index in linked list.
			/// @param Sorted Result should be sorted or not.
			void Resize(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool Sorted = true);
684 685 686 687 688 689 690 691 692
			/// Resize linked list for new interval with non-local mapping.
			/// \warning
			/// All contents of linked list will be lost after resize.
			/// @param interval_begin First index in linked list.
			/// @param interval_end Last index in linked list.
			/// @param Pre Nonlocal indices before First index in linked list.
			/// @param Post Nonlocal indices after Last index in linked list.
			/// @param Sorted Result should be sorted or not.
			void Resize(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, const std::vector<INMOST_DATA_ENUM_TYPE> & Pre, const std::vector<INMOST_DATA_ENUM_TYPE> & Post, bool Sorted = true);
Kirill Terekhov's avatar
Kirill Terekhov committed
693
#if defined(USE_SOLVER)
694 695 696 697 698
			/// Constructor that gets sizes from the matrix, including non-local mapping.
			/// @param A Matrix to get sizes from.
			/// @param Sorted Result should be sorted.
			RowMerger(const Matrix & A, bool Sorted = true);
			/// Resize linked list for new matrix, including non-local mapping.
699 700 701 702 703
			/// \warning
			/// All contents of linked list will be lost after resize.
			/// @param A Matrix to get sizes from.
			/// @param Sorted Result should be sorted or not.
			void Resize(const Matrix & A, bool Sorted = true);
Kirill Terekhov's avatar
Kirill Terekhov committed
704
#endif //USE_SOLVER
705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720
			/// Clear linked list.
			void Clear();
			/// Add a row with a coefficient into empty linked list.
			/// This routine should be a bit faster then RowMerger::AddRow
			/// for empty linked list. It may result in an unexpected behavior
			/// for non-empty linked list, asserts will fire in debug mode.
			/// @param coef Coefficient to multiply row values.
			/// @param r A row to be added.
			/// @param PreSortRow Sort values of the row before adding. Will be activated only for sorted linked lists.
			void PushRow(INMOST_DATA_REAL_TYPE coef, Row & r, bool PreSortRow);
			/// Add a row with a coefficient into non-empty linked list.
			/// Use RowMerger::PushRow for empty linked list.
			/// @param coef Coefficient to multiply row values.
			/// @param r A row to be added.
			/// @param PreSortRow Sort values of the row before adding. Will be activated only for sorted linked lists.
			void AddRow(INMOST_DATA_REAL_TYPE coef, Row & r, bool PreSortRow);
721 722 723 724 725 726 727 728 729 730 731 732 733
			/// Add a row with a coefficient into empty linked list.
			/// This routine should be a bit faster then RowMerger::AddRow
			/// for empty linked list. It may result in an unexpected behavior
			/// for non-empty linked list, asserts will fire in debug mode.
			/// @param coef Coefficient to multiply row values.
			/// @param r A row to be added.
			void PushRow(INMOST_DATA_REAL_TYPE coef, const Row & r);
			/// Add a row with a coefficient into non-empty linked list.
			/// Use RowMerger::PushRow for empty linked list.
			/// @param coef Coefficient to multiply row values.
			/// @param r A row to be added.
			void AddRow(INMOST_DATA_REAL_TYPE coef, const Row & r);
			/// Multiply all entries of linked list by a coefficient.
734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768
			/// @param coef A coefficient for multiplication.
			void Multiply(INMOST_DATA_REAL_TYPE coef);
			/// Place entries from linked list into row.
			/// \warning
			/// All contents of the row will be overwritten.
			/// If you want contents of the row to be added
			/// use AddRow with this row in advance.
			/// @param r A row to be filled.
			void RetriveRow(Row & r);
			//INMOST_DATA_REAL_TYPE ScalarProd(RowMerger & other);
			/// Get current number of nonzeros from linked list.
			INMOST_DATA_ENUM_TYPE Size() {return Nonzeros;}
			/// Retrive/add an entry from/to linked list.
			/// @param pos Position in the list.
			INMOST_DATA_REAL_TYPE & operator [] (INMOST_DATA_ENUM_TYPE pos);
			/// Retrive an entry from linked list.
			/// \warning
			/// Will fire an exception if there is no entry.
			/// @param pos Position in the list.
			INMOST_DATA_REAL_TYPE operator [] (INMOST_DATA_ENUM_TYPE pos) const;
			/// Operation of the form c = alpha a + beta b
			/// \warning
			/// Linked list must be clear before operation.
			/// @param c Row c. This will be overwritten.
			/// @param alpha Multiplier for row a.
			/// @param a Row a.
			/// @param beta Multiplier for row b.
			/// @param b Row b.
			void Merge(Row & c, INMOST_DATA_REAL_TYPE alpha, const Row & a, INMOST_DATA_REAL_TYPE beta, const Row & b)
			{
				PushRow(alpha,a);
				AddRow(beta,b);
				RetriveRow(c);
				Clear();
			}
769
			///Retrive iterator for the first element.
770
			iterator Begin() {return iterator(&LinkedList);}
771
			///Retrive iterator for the position beyond the last element.
772 773
			iterator End() {iterator ret(&LinkedList); ret.pos = EOL; return ret;}
		};
Kirill Terekhov's avatar
Kirill Terekhov committed
774
#endif //defined(USE_SOLVER) || defined(USE_AUTODIFF)
775
	} //namespace Sparse
Kirill Terekhov's avatar
Kirill Terekhov committed
776
} //namespace INMOST
Kirill Terekhov's avatar
Kirill Terekhov committed
777

778
#endif //INMOST_SPARSE_INCLUDED