inmost_solver.h 22.7 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1
2
3
4
5
6
7
8
9
10
11
#pragma once
#ifndef INMOST_SOLVER_INCLUDED
#define INMOST_SOLVER_INCLUDED


#include "inmost_common.h"
//#include "solver_prototypes.hpp"

#if defined(USE_SOLVER)
namespace INMOST
{
12
13
14
15
16
17
18
19
	/// Main class to set and solve linear system.
	/// Solver class is used to set the coefficient Matrix, the right-hand side Vector
	/// and the initial guess Vector, construct the preconditioner and Solve
	/// the linear system.
	/// Formally, Solver class is independent of INMOST::Mesh class.
	/// @see Solver::Matrix
	/// @see Solver::Vector
	/// @see Solver::Solve
Kirill Terekhov's avatar
Kirill Terekhov committed
20
21
22
23
24
	class Solver
	{
	private:
		static INMOST_MPI_Type RowEntryType; //prepared in Initialize
	public:
25
		/// Type of the Solver can be currently used in this version of INMOST.
Kirill Terekhov's avatar
Kirill Terekhov committed
26
27
		enum Type
		{
28
29
30
31
			INNER_ILU2,   ///< inner Solver based on second order ILU factorization.
			INNER_MLILUC, ///< inner Solver based on Saad multilevel ILU with pivoting.
			PETSC,        ///< external Solver PETSc, @see http://www.mcs.anl.gov/petsc/.
			ANI           ///< external Solver from ANI3D based on ILU2 (sequential Fortran version).
Kirill Terekhov's avatar
Kirill Terekhov committed
32
33
34
35
36
37
38
39
40
		};

		static INMOST_MPI_Type & GetRowEntryType() {return RowEntryType;}
		
		//solver.cpp::::::::::::::::::::::::::::::::::::::::::::::::::::
	public:
		class Matrix;
		class Vector;

41
		/// Base class for low level operations with objects of Solver class.
Kirill Terekhov's avatar
Kirill Terekhov committed
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
		class OrderInfo
		{
		private:
			typedef std::vector<INMOST_DATA_ENUM_TYPE> storage_type;
			storage_type global_to_proc; //stores ends of all non-overlapping intervals of elements, owned by this processor
			storage_type global_overlap; //stores pairs: [begin,end) of overlapping intervals of rows
			std::vector<INMOST_DATA_ENUM_TYPE> vector_exchange_recv, vector_exchange_send;
			std::vector<INMOST_DATA_REAL_TYPE> send_storage, recv_storage;
			std::vector<INMOST_MPI_Request> send_requests, recv_requests;
			std::vector<INMOST_DATA_ENUM_TYPE> extended_indexes;

			//remote indexes
			INMOST_DATA_ENUM_TYPE local_vector_begin, local_vector_end;
			INMOST_DATA_ENUM_TYPE initial_matrix_begin, initial_matrix_end; //local interval of matrix
			INMOST_DATA_ENUM_TYPE local_matrix_begin, local_matrix_end; //local interval of matrix

			bool have_matrix;
			INMOST_MPI_Comm comm;
			int rank,size;
		public:
Kirill Terekhov's avatar
Kirill Terekhov committed
62
			void Clear();
63
			/// Return true if Matrix data have already been specified.
Kirill Terekhov's avatar
Kirill Terekhov committed
64
65
66
67
68
			bool & HaveMatrix() { return have_matrix; }
			OrderInfo();
			OrderInfo(const OrderInfo & other);
			OrderInfo & operator =(OrderInfo const & other);
			~OrderInfo();
69
70
71
72
73
			/// Prepare parallel state of the Matrix with specified overlap size.
			/// This state of the matrix can be used, for instance, to construct
			/// the preconditioner for Additive Swartz method.
			/// @param m Matrix to be expanded.
			/// @param overlap Overlap size, viz. the number of overlap layers.
Kirill Terekhov's avatar
Kirill Terekhov committed
74
			void PrepareMatrix(Matrix & m, INMOST_DATA_ENUM_TYPE overlap);
75
			/// Restore initial nonparallel state of the Matrix with no overlap.
Kirill Terekhov's avatar
Kirill Terekhov committed
76
			void RestoreMatrix(Matrix & m);
77
			/// Prepare parallel state of the Vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
78
			void PrepareVector(Vector & v);
79
			/// Restore initial nonparallel state of the Vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
80
			void RestoreVector(Vector & v); //Restore initial nonparallel state of the vector
81
			/// Retrieve the processor number by binary search for the specified global index.
Kirill Terekhov's avatar
Kirill Terekhov committed
82
83
			INMOST_DATA_ENUM_TYPE GetProcessor(INMOST_DATA_ENUM_TYPE gind) const; //retrive processor by binary search in global_to_proc
			void      GetOverlapRegion(INMOST_DATA_ENUM_TYPE proc, INMOST_DATA_ENUM_TYPE & mbeg, INMOST_DATA_ENUM_TYPE & mend) const;
84
			/// Get the local index region for the specified process.
Kirill Terekhov's avatar
Kirill Terekhov committed
85
			void        GetLocalRegion(INMOST_DATA_ENUM_TYPE proc, INMOST_DATA_ENUM_TYPE & mbeg, INMOST_DATA_ENUM_TYPE & mend) const;
86
			/// Get the local index region for the current process.
Kirill Terekhov's avatar
Kirill Terekhov committed
87
			void       GetVectorRegion(INMOST_DATA_ENUM_TYPE & mbeg, INMOST_DATA_ENUM_TYPE & mend) const {mbeg = local_vector_begin; mend = local_vector_end;}
88
			/// Get the rank of the current communicator, i.e. the current process index.
Kirill Terekhov's avatar
Kirill Terekhov committed
89
			INMOST_DATA_ENUM_TYPE GetRank() const {return rank;}
90
			/// Get the size of the current communicator, i.e. the total number of processes used.
Kirill Terekhov's avatar
Kirill Terekhov committed
91
			INMOST_DATA_ENUM_TYPE GetSize() const {return size;}
92
			/// Update the shared data in parallel vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
93
			void                  Update    (Vector & x); // update parallel vector
94
			/// Sum shared values in parallel vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
95
			void                  Accumulate(Vector & x); // sum shared values in parallel vector
96
			/// Get the sum of num elements of real array on all processes.
Kirill Terekhov's avatar
Kirill Terekhov committed
97
			void                  Integrate(INMOST_DATA_REAL_TYPE * inout, INMOST_DATA_ENUM_TYPE num);
98
			/// Get the communicator which the solver is associated with.
Kirill Terekhov's avatar
Kirill Terekhov committed
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
			INMOST_MPI_Comm         GetComm() {return comm;}
			// Access to arrays below allows to organize manual exchange
			INMOST_MPI_Request    * GetSendRequests() {assert(!send_requests.empty()); return &send_requests[0];}
			INMOST_MPI_Request    * GetRecvRequests() {assert(!recv_requests.empty()); return &recv_requests[0];}
			INMOST_DATA_ENUM_TYPE   GetSendRequestsSize() {return static_cast<INMOST_DATA_ENUM_TYPE>(send_requests.size());}
			INMOST_DATA_ENUM_TYPE   GetRecvRequestsSize() {return static_cast<INMOST_DATA_ENUM_TYPE>(recv_requests.size());}
			INMOST_DATA_ENUM_TYPE * GetSendExchangeArray() {assert(!vector_exchange_send.empty()); return &vector_exchange_send[0];}
			INMOST_DATA_ENUM_TYPE   GetSendExchangeSize() {return static_cast<INMOST_DATA_ENUM_TYPE>(send_storage.size());}
			INMOST_DATA_ENUM_TYPE * GetRecvExchangeArray() {assert(!vector_exchange_recv.empty()); return &vector_exchange_recv[0];}
			INMOST_DATA_ENUM_TYPE   GetRecvExchangeSize() {return static_cast<INMOST_DATA_ENUM_TYPE>(recv_storage.size());}
			//for debug
			//~ void                 BeginSequentialCode() {for(int i = 0; i < rank; i++) MPI_Barrier(comm);}
			//~ void                   EndSequentialCode() {for(int i = rank; i < size; i++) MPI_Barrier(comm);}
		};

114
115
116
117
		/// 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 Solver::Solve
Kirill Terekhov's avatar
Kirill Terekhov committed
118
119
120
121
122
123
124
125
126
127
128
129
		class Vector
		{
		public:
			typedef interval<INMOST_DATA_ENUM_TYPE,INMOST_DATA_REAL_TYPE> Entries;
			typedef Entries::iterator iterator;
			typedef Entries::const_iterator const_iterator;
		private:
			INMOST_MPI_Comm comm;
			Entries data;
			std::string name;
			bool is_parallel;
		public:
130
131
132
133
134
			/// 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.
Kirill Terekhov's avatar
Kirill Terekhov committed
135
136
137
138
			Vector(std::string _name = "", INMOST_DATA_ENUM_TYPE start = 0, INMOST_DATA_ENUM_TYPE end = 0, INMOST_MPI_Comm _comm = INMOST_MPI_COMM_WORLD);
			Vector(const Vector & other);
			Vector & operator =(Vector const & other);
			~Vector();
139
			/// Return reference to i-th element of the vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
140
			INMOST_DATA_REAL_TYPE & operator [](INMOST_DATA_ENUM_TYPE i) {return data[i];}
141
			/// Return i-th element of the vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
142
			INMOST_DATA_REAL_TYPE   operator [](INMOST_DATA_ENUM_TYPE i) const {return data[i];}
143
			/// Return the global size of the vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
144
145
146
147
148
149
			INMOST_DATA_ENUM_TYPE  Size() const { return static_cast<INMOST_DATA_ENUM_TYPE>(data.size()); }
			iterator             Begin() {return data.begin();}
			const_iterator       Begin() const {return data.begin();}
			iterator             End() {return data.end();}
			const_iterator       End() const {return data.end();}
			bool                 Empty() const {return data.empty();}
150
			/// Set the start and the end of the distributed vector interval.
Kirill Terekhov's avatar
Kirill Terekhov committed
151
			void                 SetInterval(INMOST_DATA_ENUM_TYPE   start, INMOST_DATA_ENUM_TYPE   end)       {data.set_interval_beg(start); data.set_interval_end(end);}
152
			/// Get the start and the end of the distributed vector interval.
Kirill Terekhov's avatar
Kirill Terekhov committed
153
154
			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);}
155
			/// Get the first index of the distributed vector interval.
Kirill Terekhov's avatar
Kirill Terekhov committed
156
			INMOST_DATA_ENUM_TYPE  GetFirstIndex() const {return data.get_interval_beg();}
157
			/// Get the communicator which the vector is associated with.
Kirill Terekhov's avatar
Kirill Terekhov committed
158
159
			INMOST_MPI_Comm        GetCommunicator() const {return comm;}

160
			/// Get the scalar product of the specified interval of the distributed vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
161
162
163
			INMOST_DATA_REAL_TYPE ScalarProd(Vector const & other, INMOST_DATA_ENUM_TYPE index_begin, INMOST_DATA_ENUM_TYPE index_end) const;


164
			/// Save the distributed vector to a single data file using parallel MPI I/O.
Kirill Terekhov's avatar
Kirill Terekhov committed
165
			void                 Save(std::string file);
166
167
168
			/// 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).
Kirill Terekhov's avatar
Kirill Terekhov committed
169
170
171
			void                 Load(std::string file, INMOST_DATA_ENUM_TYPE mbeg = ENUMUNDEF, INMOST_DATA_ENUM_TYPE mend = ENUMUNDEF);
			
			bool                 & isParallel() {return is_parallel;}
172
			/// Get the vector name specified in the main constructor.
Kirill Terekhov's avatar
Kirill Terekhov committed
173
			std::string          GetName() {return name;}
Kirill Terekhov's avatar
Kirill Terekhov committed
174

175
			/// Clear all data of the current vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
176
			void Clear() {data.clear();}
Kirill Terekhov's avatar
Kirill Terekhov committed
177
178
179
180
181
			//~ friend class Solver;
		};



182
		/// Class to store the sparse matrix row.
Kirill Terekhov's avatar
Kirill Terekhov committed
183
184
185
186
		class Row 
		{
		public:
			
187
			/// Entry of the sparse matrix row.
Kirill Terekhov's avatar
Kirill Terekhov committed
188
189
			typedef struct entry_s 
			{
190
191
				INMOST_DATA_ENUM_TYPE first;  ///< the column number of the row element.
				INMOST_DATA_REAL_TYPE second; ///< the real value of the row element.
Kirill Terekhov's avatar
Kirill Terekhov committed
192
				entry_s() :first(0), second(0.0) {}
Kirill Terekhov's avatar
Kirill Terekhov committed
193
				entry_s(const entry_s & other) :first(other.first), second(other.second) {}//{std::cout << __FUNCTION__ << " " << first << " " << second << std::endl;}
Kirill Terekhov's avatar
Kirill Terekhov committed
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
				entry_s(INMOST_DATA_ENUM_TYPE first, INMOST_DATA_REAL_TYPE second):first(first),second(second){}
				entry_s & operator =(entry_s const & other) {first = other.first, second = other.second; return *this;}
				bool operator < (const entry_s & other) const { return first < other.first || (first == other.first && second < other.second); }
			} entry;
			
			typedef dynarray<entry,16> Entries; //replace later with more memory-efficient chunk_array, with first chunk in stack
			//typedef array<entry> Entries;
			//typedef std::vector<entry> Entries;
			//typedef sparse_data<INMOST_DATA_ENUM_TYPE,INMOST_DATA_REAL_TYPE> Entries;
			//typedef Entries::pair entry; //for sparse_data
			typedef Entries::iterator iterator;
			typedef Entries::const_iterator const_iterator;
			typedef Entries::reverse_iterator reverse_iterator;
			typedef Entries::const_reverse_iterator const_reverse_iterator;

			bool modified_pattern; //remove this in future
		private:
#if defined(USE_OMP)
			omp_lock_t lock;
#endif
			bool marker;
			Entries data;
		public:
Kirill Terekhov's avatar
Kirill Terekhov committed
217
			void Report() {data.report_addr();}
Kirill Terekhov's avatar
Kirill Terekhov committed
218
219
220
221
222
223
224
225
226
227
228
229
			void SetMarker() { marker = true; }
			void RemMarker() { marker = false; }
			bool GetMarker() { return marker; }
			Row() :data() 
			{
#if defined(USE_OMP)
				omp_init_lock(&lock);
#endif
				modified_pattern = marker = false;
			}
			Row(const Row & other) :marker(other.marker),data(other.data) 
			{ 
Kirill Terekhov's avatar
Kirill Terekhov committed
230
231
232
				//std::cout << __FUNCTION__ << " ";
				//for(iterator it = Begin(); it != End(); ++it) std::cout << it->first << "," << it->second << " ";
				//std::cout << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
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
#if defined(USE_OMP)
				omp_init_lock(&lock);
#endif
				modified_pattern = other.modified_pattern; 
			}
			Row(entry * pbegin, entry * pend) :data(pbegin, pend) 
			{ 
#if defined(USE_OMP)
				omp_init_lock(&lock);
#endif
				modified_pattern = true; marker = false; 
			}
			void Lock() 
			{
#if defined(USE_OMP)
				omp_set_lock(&lock); 
#endif
			}
			void Unlock() 
			{ 
#if defined(USE_OMP)
				omp_unset_lock(&lock); 
#endif
			}
			~Row() {}
			Row & operator = (Row const & other) { data = other.data; marker = other.marker; return *this; }
259
			/// The operator [] used to fill the sparse matrix row, but not to access individual elements of the row.
Kirill Terekhov's avatar
Kirill Terekhov committed
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
			INMOST_DATA_REAL_TYPE & operator [](INMOST_DATA_ENUM_TYPE i) // use to fill matrix, not to access individual elements
			{
				//for sparse_data type
				//return data[i];
				//for dynarray or array
				
				for(unsigned it = 0; it < data.size(); ++it)
					if( data[it].first == i ) return data[it].second;
				entry new_entry;
				new_entry.first = i;
				new_entry.second = 0;
				data.push_back(new_entry);
				modified_pattern = true;
				return data.back().second;
				
			}
276
			/// The operator [] used to access individual elements of the row.
Kirill Terekhov's avatar
Kirill Terekhov committed
277
278
279
280
281
282
283
284
285
286
287
288
			INMOST_DATA_REAL_TYPE operator[](INMOST_DATA_ENUM_TYPE i) const
			{
				//for sparse data type
				//return data[i];

				for (unsigned 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;
			}
			//void           Reserve(INMOST_DATA_ENUM_TYPE num) { data.reserve(num);}
289
			/// Clear all data of the current row.
Kirill Terekhov's avatar
Kirill Terekhov committed
290
291
			void                  Clear() { data.clear(); }
			void                  Swap(Solver::Row & other) { data.swap(other.data); bool tmp = marker; marker = other.marker; other.marker = tmp; }
292
			/// The size of the sparse row, i.e. the total number of nonzero elements.
Kirill Terekhov's avatar
Kirill Terekhov committed
293
294
295
296
297
298
299
300
301
302
303
304
305
306
			INMOST_DATA_ENUM_TYPE   Size() const { return static_cast<INMOST_DATA_ENUM_TYPE>(data.size()); }
			INMOST_DATA_ENUM_TYPE & 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;}
			INMOST_DATA_ENUM_TYPE   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(); }
307
			/// Return the scalar product of the current sparse row by a dense Vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
308
309
			INMOST_DATA_REAL_TYPE   RowVec(Vector & x) const; // returns A(row) * x
			void                  MoveRow(Row & new_pos) {data = new_pos.data;} //here move constructor and std::move may be used in future
310
			/// Set the vector entries by zeroes.
Kirill Terekhov's avatar
Kirill Terekhov committed
311
			void                  Zero() {for(iterator it = Begin(); it != End(); ++it) it->second = 0;}
312
313
314
			/// Push specified element into sparse row.
			/// This function should be used only if the index is not repeated in the row
			/// and all indexes are inserted in increasing order.
Kirill Terekhov's avatar
Kirill Terekhov committed
315
316
317
			void                  Push(INMOST_DATA_ENUM_TYPE ind, INMOST_DATA_REAL_TYPE val) {data.push_back(entry(ind,val));}
		};
		
318
319
320
		/// 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
Kirill Terekhov's avatar
Kirill Terekhov committed
321
322
323
324
325
326
327
328
329
330
331
332
		class Matrix
		{
		public:
			typedef interval<INMOST_DATA_ENUM_TYPE,Solver::Row> Rows;
			typedef Rows::iterator iterator;
			typedef Rows::const_iterator const_iterator;
		private:
			INMOST_MPI_Comm comm;
			Rows data;
			std::string name;
			bool is_parallel;
		public:
333
334
335
336
337
			/// 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.
Kirill Terekhov's avatar
Kirill Terekhov committed
338
339
340
341
342
			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();
			
343
			/// Return reference to i-th Row of the matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
344
			Row & operator [](INMOST_DATA_ENUM_TYPE i) {return data[i];}
345
			/// Return reference to i-th Row of the matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
346
			const Row & operator [](INMOST_DATA_ENUM_TYPE i) const {return data[i];}
347
			/// Return the total number of rows in the matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
348
349
350
351
352
353
			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();}
354
			/// Set the start and the end row numbers of the distributed matrix interval.
Kirill Terekhov's avatar
Kirill Terekhov committed
355
			void                 SetInterval(INMOST_DATA_ENUM_TYPE   start, INMOST_DATA_ENUM_TYPE   end)       {data.set_interval_beg(start); data.set_interval_end(end);}
356
			/// Get the start and the end row numbers of the distributed matrix interval.
Kirill Terekhov's avatar
Kirill Terekhov committed
357
358
			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);}
359
			/// Get the first row index of the distributed matrix interval.
Kirill Terekhov's avatar
Kirill Terekhov committed
360
			INMOST_DATA_ENUM_TYPE  GetFirstIndex() const {return data.get_interval_beg();}
361
			/// Get the communicator which the matrix is associated with.
Kirill Terekhov's avatar
Kirill Terekhov committed
362
363
364
365
366
367
368
369
370
371
372
373
374
375
			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(Solver::Matrix & other) 
			{
				data.swap(other.data);
				name.swap(other.name);
				INMOST_MPI_Comm ctmp = comm;
				comm = other.comm;
				other.comm = ctmp;
				bool ptmp = is_parallel;
				is_parallel = other.is_parallel;
				other.is_parallel = ptmp;
			}
			
376
377
378
379
			/// Matrix-vector product of the form: y = alpha*A*x + beta * y.
			/// @param y Input/output vector.
			/// @see Solver::Vector::Zero
			void MatVec(INMOST_DATA_REAL_TYPE alpha, Solver::Vector & x, INMOST_DATA_REAL_TYPE beta, Solver::Vector & y) const; //y = alpha*A*x + beta * y
Kirill Terekhov's avatar
Kirill Terekhov committed
380

381
			/// Clear all data of the matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
382
			void Clear() {for(Matrix::iterator it = Begin(); it != End(); ++it) it->Clear(); data.clear();}
Kirill Terekhov's avatar
Kirill Terekhov committed
383

384
385
386
			/// 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).
Kirill Terekhov's avatar
Kirill Terekhov committed
387
			void				 Load(std::string file, INMOST_DATA_ENUM_TYPE beg = ENUMUNDEF, INMOST_DATA_ENUM_TYPE end = ENUMUNDEF);
388
389
			/// 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
Kirill Terekhov's avatar
Kirill Terekhov committed
390
391
392
			void                 Save(std::string file);
			bool &               isParallel() { return is_parallel; }
			
393
			/// Get the matrix name specified in the main constructor.
Kirill Terekhov's avatar
Kirill Terekhov committed
394
395
396
397
398
			std::string          GetName() {return name;}
			//~ friend class Solver;
		};
		
	private:
Kirill Terekhov's avatar
Kirill Terekhov committed
399
		static bool is_initialized, is_finalized;
Kirill Terekhov's avatar
Kirill Terekhov committed
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
		INMOST_MPI_Comm comm;
		std::string name;
		INMOST_DATA_ENUM_TYPE local_size, global_size;
		INMOST_DATA_ENUM_TYPE last_it;
		INMOST_DATA_REAL_TYPE last_resid;
		OrderInfo info;
		INMOST_DATA_ENUM_TYPE overlap;
		
		void * solver_data;
		void * matrix_data;
		void * precond_data;
		
		void * rhs_data;
		void * solution_data;
		
		Type _pack;
		Solver(const Solver & other);// prohibit copy
		Solver & operator =(Solver const & other); //prohibit assignment
	public:
419
		/// Set the solver parameter of the integer type.
Kirill Terekhov's avatar
Kirill Terekhov committed
420
		void SetParameterEnum(std::string name, INMOST_DATA_ENUM_TYPE value);
421
		/// Set the solver parameter of the real type.
Kirill Terekhov's avatar
Kirill Terekhov committed
422
		void SetParameterReal(std::string name, INMOST_DATA_REAL_TYPE value);
423
		/// Get the used defined name of the Solver.
Kirill Terekhov's avatar
Kirill Terekhov committed
424
425
		std::string          GetName() {return name;}
	
426
		/// Get the package Type.
Kirill Terekhov's avatar
Kirill Terekhov committed
427
428
		Type GetPackage() const {return _pack;}
		
429
430
		/// Return the number of iterations performed by the last solution.
		/// @see Solver::Solve
Kirill Terekhov's avatar
Kirill Terekhov committed
431
		INMOST_DATA_ENUM_TYPE Iterations();
432
433
		/// Return the final residual achieved by the last solution.
		/// @see Solver::Solve
Kirill Terekhov's avatar
Kirill Terekhov committed
434
		INMOST_DATA_REAL_TYPE Residual();
435
436
437
438
		/// Set the matrix and construct the preconditioner.
		/// @param OldPreconditioner If this parameter is set to true,
		/// then the previous preconditioner will be used,
		/// otherwise the new preconditioner will be constructed. 
Kirill Terekhov's avatar
Kirill Terekhov committed
439
		void SetMatrix(Matrix & A, bool OldPreconditioner = false);
440
441
442
443
444
445
		/// Solver the linear system: A*x = b.
		/// @param RHS The right-hand side Vector b.
		/// @param SOL The initial guess to the solution on input and the solution Vector x on return.
		/// It is assumed that the coefficient matrix A have been set
		/// and the preconditioner have been already constructed.
		/// @see Solver::SetMatrix
Kirill Terekhov's avatar
Kirill Terekhov committed
446
		bool Solve(Vector & RHS, Vector & SOL);
447
448
		/// Get the reason of convergence or divergence of the last solution.
		/// @see Solver::Solve
Kirill Terekhov's avatar
Kirill Terekhov committed
449
		std::string GetReason();
Kirill Terekhov's avatar
Kirill Terekhov committed
450

451
452
453
454
455
456
457
458
		/// Main constructor of the solver.
		/// @param pack The package Type to be used for solution.
		/// @param _name The user specified name of the current solver.
		/// @param _comm Communicator for parallel data exchanges, MPI_COMM_WORLD by default.
		/// @see Solver::Initialize
		/// @see Solver::SetMatrix
		/// @see Solver::Solve
		/// @see Solver::Finalize
Kirill Terekhov's avatar
Kirill Terekhov committed
459
460
		Solver(Type pack, std::string _name = "", INMOST_MPI_Comm comm = INMOST_MPI_COMM_WORLD);
		~Solver();
461
462
463
464
465
466
467
468
		/// Initialize the stage of parallel solution.
		/// If MPI is not initialized yet, then it will be initialized.
		/// @param argc The number of arguments transmitted to the function main.
		/// @param argv The pointer to arguments transmitted to the function main.
		/// @param database Usually the name of the file with the Solver parameters.
		/// The shortest call to this function with the default solver parameters is the following: Initialize(NULL,NULL,"");
		/// @see Solver::Finalize
		/// @see Solver::isInitialized
Kirill Terekhov's avatar
Kirill Terekhov committed
469
		static void Initialize(int * argc, char *** argv, const char * database = "");
470
471
472
473
474
		/// Finalize the stage of parallel solution.
		/// If MPI was initialized in Solver::Initialize, then it will be finalized.
		/// By this reason, do not use any MPI function after call to this function.
		/// @see Solver::Initialize
		/// @see Solver::isFinalized
Kirill Terekhov's avatar
Kirill Terekhov committed
475
		static void Finalize();
Kirill Terekhov's avatar
Kirill Terekhov committed
476
477
478
		static bool isInitialized() {return is_initialized;}
		static bool isFinalized() {return is_finalized;}

479
		/// Clear all internal data of the current solver including matrix, preconditioner etc.
Kirill Terekhov's avatar
Kirill Terekhov committed
480
		void Clear();
Kirill Terekhov's avatar
Kirill Terekhov committed
481
482
483
484
485
486
	};
}

#endif // USE_SOLVER

#endif // INMOST_SOLVER_INCLUDED