inmost_solver.h 37.7 KB
Newer Older
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
1

Kirill Terekhov's avatar
Kirill Terekhov committed
2
3
4
5
6
7
8
#ifndef INMOST_SOLVER_INCLUDED
#define INMOST_SOLVER_INCLUDED


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

Kirill Terekhov's avatar
Kirill Terekhov committed
9
#define DEFAULT_ADDITIVE_SCHWARTZ_OVERLAP             1
Kirill Terekhov's avatar
Kirill Terekhov committed
10
11
#define DEFAULT_ABSOLUTE_TOLERANCE                    1.0e-5
#define DEFAULT_RELATIVE_TOLERANCE                    1.0e-12
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
12
#define DEFAULT_DIVERGENCE_TOLERANCE                  1.0e+100
Kirill Terekhov's avatar
Kirill Terekhov committed
13
#define DEFAULT_MAXIMUM_ITERATIONS                    2500
Kirill Terekhov's avatar
Kirill Terekhov committed
14
15
16
#define DEFAULT_SOLVER_GMRES_SUBSTEPS                 2
#define DEFAULT_PRECONDITIONER_DROP_TOLERANCE         0.005
#define DEFAULT_PRECONDITIONER_REUSE_TOLERANCE        0.00005
Kirill Terekhov's avatar
Kirill Terekhov committed
17
#define DEFAULT_PRECONDITIONER_FILL_LEVEL             3
Kirill Terekhov's avatar
Kirill Terekhov committed
18
19
20
#define DEFAULT_PRECONDITIONER_DDPQ_TOLERANCE         0.75
#define DEFAULT_PRECONDITIONER_REORDER_NONZEROS       1
#define DEFAULT_PRECONDITIONER_RESCALE_ITERS          6
Kirill Terekhov's avatar
Kirill Terekhov committed
21
22
23
#define DEFAULT_PRECONDITIONER_CONDITION_ESTIMATION   1
#define DEFAULT_PRECONDITIONER_ADAPT_DDPQ_TOLERANCE   1

Kirill Terekhov's avatar
Kirill Terekhov committed
24
25
26
#if defined(USE_SOLVER)
namespace INMOST
{
27
28
29
30
	/// 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.
31
	///
32
33
34
35
	/// 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
36
37
38
39
40
	class Solver
	{
	private:
		static INMOST_MPI_Type RowEntryType; //prepared in Initialize
	public:
41
		/// Type of the Solver can be currently used in this version of INMOST.
Kirill Terekhov's avatar
Kirill Terekhov committed
42
43
		enum Type
		{
44
			INNER_ILU2,     ///< inner Solver based on BiCGStab(L) solver with second order ILU factorization as preconditioner.
45
46
47
48
49
50
51
52
53
			INNER_DDPQILUC, ///< inner Solver based on BiCGStab(L) solver with second order Crout-ILU with inversed-based condition estimation and unsymmetric reordering for diagonal dominance as preconditioner.
			INNER_MPTILUC,  ///< inner Solver based on BiCGStab(L) solver with second order Crout-ILU with inversed-based condition estimation and maximum product transversal reordering as preconditioner.
			INNER_MPTILU2,  ///< inner Solver based on BiCGStab(L) solver with second order ILU and maximum product transversal reordering as preconditioner.
			Trilinos_Aztec, ///< external Solver AztecOO from Trilinos package.
			Trilinos_Belos, ///< external Solver Belos from Trilinos package, currently without preconditioner.
			Trilinos_ML,    ///< external Solver AztecOO with ML preconditioner.
			Trilinos_Ifpack,///< external Solver AztecOO with Ifpack preconditioner.
			PETSc,          ///< external Solver PETSc, @see http://www.mcs.anl.gov/petsc/
			ANI,            ///< external Solver from ANI3D based on ILU2 (sequential Fortran version), @see http://ani3d.sourceforge.net/
54
55
			FCBIILU2,       ///< external FCBIILU2 Solver (BIILU2 parallel F2C version).
			K3BIILU2        ///< inner    K3BIILU2 Solver (BIILU2 parallel version).
Kirill Terekhov's avatar
Kirill Terekhov committed
56
57
58
59
60
61
62
63
64
		};

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

65
		/// Base class for low level operations with objects of Solver class.
Kirill Terekhov's avatar
Kirill Terekhov committed
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
		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
86
			void Clear();
87
			/// Return true if Matrix data have already been specified.
Kirill Terekhov's avatar
Kirill Terekhov committed
88
89
90
91
92
			bool & HaveMatrix() { return have_matrix; }
			OrderInfo();
			OrderInfo(const OrderInfo & other);
			OrderInfo & operator =(OrderInfo const & other);
			~OrderInfo();
93
94
95
96
97
			/// 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
98
			void PrepareMatrix(Matrix & m, INMOST_DATA_ENUM_TYPE overlap);
99
			/// Restore initial nonparallel state of the Matrix with no overlap.
Kirill Terekhov's avatar
Kirill Terekhov committed
100
			void RestoreMatrix(Matrix & m);
101
			/// Prepare parallel state of the Vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
102
			void PrepareVector(Vector & v) const;
103
			/// Restore initial nonparallel state of the Vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
104
			void RestoreVector(Vector & v) const;
105
			/// Retrieve the processor number by binary search for the specified global index.
Igor Konshin's avatar
Igor Konshin committed
106
			INMOST_DATA_ENUM_TYPE GetProcessor(INMOST_DATA_ENUM_TYPE gind) const; //retrieve processor by binary search in global_to_proc
Kirill Terekhov's avatar
Kirill Terekhov committed
107
			void      GetOverlapRegion(INMOST_DATA_ENUM_TYPE proc, INMOST_DATA_ENUM_TYPE & mbeg, INMOST_DATA_ENUM_TYPE & mend) const;
108
			/// Get the local index region for the specified process.
Kirill Terekhov's avatar
Kirill Terekhov committed
109
			void        GetLocalRegion(INMOST_DATA_ENUM_TYPE proc, INMOST_DATA_ENUM_TYPE & mbeg, INMOST_DATA_ENUM_TYPE & mend) const;
110
			/// Get the local index region for the current process.
Kirill Terekhov's avatar
Kirill Terekhov committed
111
			void       GetVectorRegion(INMOST_DATA_ENUM_TYPE & mbeg, INMOST_DATA_ENUM_TYPE & mend) const {mbeg = local_vector_begin; mend = local_vector_end;}
112
			/// Get the rank of the current communicator, i.e. the current process index.
Kirill Terekhov's avatar
Kirill Terekhov committed
113
			INMOST_DATA_ENUM_TYPE GetRank() const {return rank;}
114
			/// Get the size of the current communicator, i.e. the total number of processes used.
Kirill Terekhov's avatar
Kirill Terekhov committed
115
			INMOST_DATA_ENUM_TYPE GetSize() const {return size;}
116
			/// Update the shared data in parallel vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
117
			void                  Update    (Vector & x); // update parallel vector
118
			/// Sum shared values in parallel vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
119
			void                  Accumulate(Vector & x); // sum shared values in parallel vector
120
			/// Get the sum of num elements of real array on all processes.
Kirill Terekhov's avatar
Kirill Terekhov committed
121
			void                  Integrate(INMOST_DATA_REAL_TYPE * inout, INMOST_DATA_ENUM_TYPE num) const;
122
			/// Get the communicator which the solver is associated with.
Kirill Terekhov's avatar
Kirill Terekhov committed
123
			INMOST_MPI_Comm         GetComm() const {return comm;}
Kirill Terekhov's avatar
Kirill Terekhov committed
124
125
126
127
128
129
130
131
132
133
134
135
			// 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);}
Kirill Terekhov's avatar
Kirill Terekhov committed
136

Kirill Terekhov's avatar
Kirill Terekhov committed
137
138
139
			// Get the scalar product of the specified interval of the distributed vector.
			// Conflicts with OpenMP, should not be used in future
			//void ScalarProd(Vector const & left, Vector const & right, INMOST_DATA_ENUM_TYPE index_begin, INMOST_DATA_ENUM_TYPE index_end, INMOST_DATA_REAL_TYPE & sum) const;
Kirill Terekhov's avatar
Kirill Terekhov committed
140
141
		};

142
143
144
145
		/// 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
146
147
148
149
150
151
152
153
154
155
156
157
		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:
158
159
160
161
162
			/// 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
163
164
165
166
			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();
167
			/// Return reference to i-th element of the vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
168
			INMOST_DATA_REAL_TYPE & operator [](INMOST_DATA_ENUM_TYPE i) {return data[i];}
169
			/// Return i-th element of the vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
170
			INMOST_DATA_REAL_TYPE   operator [](INMOST_DATA_ENUM_TYPE i) const {return data[i];}
171
			/// Return the global size of the vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
172
173
174
175
176
177
			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();}
178
			/// Set the start and the end of the distributed vector interval.
Kirill Terekhov's avatar
Kirill Terekhov committed
179
			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);}
180
			/// Get the start and the end of the distributed vector interval.
Kirill Terekhov's avatar
Kirill Terekhov committed
181
182
			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);}
183
			/// Get the first index of the distributed vector interval.
Kirill Terekhov's avatar
Kirill Terekhov committed
184
			INMOST_DATA_ENUM_TYPE  GetFirstIndex() const {return data.get_interval_beg();}
185
			/// Get the communicator which the vector is associated with.
Kirill Terekhov's avatar
Kirill Terekhov committed
186
187
			INMOST_MPI_Comm        GetCommunicator() const {return comm;}

Kirill Terekhov's avatar
Kirill Terekhov committed
188
			void Swap(Solver::Vector & other) {data.swap(other.data); name.swap(other.name); std::swap(is_parallel,other.is_parallel); std::swap(comm,other.comm);}
Kirill Terekhov's avatar
Kirill Terekhov committed
189
190


191
			/// Save the distributed vector to a single data file using parallel MPI I/O.
Kirill Terekhov's avatar
Kirill Terekhov committed
192
			void                 Save(std::string file);
193
194
195
			/// 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
196
197
198
			void                 Load(std::string file, INMOST_DATA_ENUM_TYPE mbeg = ENUMUNDEF, INMOST_DATA_ENUM_TYPE mend = ENUMUNDEF);
			
			bool                 & isParallel() {return is_parallel;}
199
			/// Get the vector name specified in the main constructor.
Kirill Terekhov's avatar
Kirill Terekhov committed
200
			std::string          GetName() {return name;}
Kirill Terekhov's avatar
Kirill Terekhov committed
201

202
			/// Clear all data of the current vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
203
			void Clear() {data.clear();}
Kirill Terekhov's avatar
Kirill Terekhov committed
204
205
206
207
208
			//~ friend class Solver;
		};



209
		/// Class to store the sparse matrix row.
Kirill Terekhov's avatar
Kirill Terekhov committed
210
211
212
213
		class Row 
		{
		public:
			
214
			/// Entry of the sparse matrix row.
Kirill Terekhov's avatar
Kirill Terekhov committed
215
216
			typedef struct entry_s 
			{
217
218
				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
Fixes    
Kirill Terekhov committed
219
220
221
222
				//entry_s() :first(0), second(0.0) {}
				//entry_s(const entry_s & other) :first(other.first), second(other.second) {}//{std::cout << __FUNCTION__ << " " << first << " " << second << std::endl;}
				//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;}
Kirill Terekhov's avatar
Kirill Terekhov committed
223
224
				bool operator < (const entry_s & other) const { return first < other.first || (first == other.first && second < other.second); }
			} entry;
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
225
226
227
228
229
230
231
			__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;
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
232
233

		private:
Kirill Terekhov's avatar
Kirill Terekhov committed
234
235
236
237
238
			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
Kirill Terekhov's avatar
Kirill Terekhov committed
239
240
			
		public:
Kirill Terekhov's avatar
Kirill Terekhov committed
241
242
243
244
245
246
247
248
249
250
251
252
253
			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
254
			void Report() {data.report_addr();}
Kirill Terekhov's avatar
Kirill Terekhov committed
255
256
257
258
259
260
261
262
263
264
265
266
			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
267
268
269
				//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
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
#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; }
296
			/// 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
297
298
299
300
301
302
			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
				
Kirill Terekhov's avatar
Kirill Terekhov committed
303
				for(Entries::size_type it = 0; it < data.size(); ++it)
Kirill Terekhov's avatar
Kirill Terekhov committed
304
305
306
307
308
309
310
311
312
					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;
				
			}
313
			/// The operator [] used to access individual elements of the row.
Kirill Terekhov's avatar
Kirill Terekhov committed
314
315
316
317
318
			INMOST_DATA_REAL_TYPE operator[](INMOST_DATA_ENUM_TYPE i) const
			{
				//for sparse data type
				//return data[i];

Kirill Terekhov's avatar
Kirill Terekhov committed
319
				for (Entries::size_type it = 0; it < data.size(); ++it) if (data[it].first == i) return data[it].second;
Kirill Terekhov's avatar
Kirill Terekhov committed
320
321
322
323
324
325

				//you should not come here
				assert(false);
				return 1.0e20;
			}
			//void           Reserve(INMOST_DATA_ENUM_TYPE num) { data.reserve(num);}
326
			/// Clear all data of the current row.
Kirill Terekhov's avatar
Kirill Terekhov committed
327
328
			void                    Clear() { data.clear(); }
			void                    Swap(Solver::Row & other) { data.swap(other.data); bool tmp = marker; marker = other.marker; other.marker = tmp; }
329
			/// The size of the sparse row, i.e. the total number of nonzero elements.
Kirill Terekhov's avatar
Kirill Terekhov committed
330
331
332
333
334
335
			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;}
			
Kirill Terekhov's avatar
Kirill Terekhov committed
336
337
338
339
340
341
342
343
			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(); }
344
			/// Return the scalar product of the current sparse row by a dense Vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
345
			INMOST_DATA_REAL_TYPE   RowVec(Vector & x) const; // returns A(row) * x
Kirill Terekhov's avatar
Kirill Terekhov committed
346
			void                    MoveRow(Row & new_pos) {data = new_pos.data;} //here move constructor and std::move may be used in future
347
			/// Set the vector entries by zeroes.
Kirill Terekhov's avatar
Kirill Terekhov committed
348
			void                    Zero() {for(iterator it = Begin(); it != End(); ++it) it->second = 0;}
349
			/// Push specified element into sparse row.
Kirill Terekhov's avatar
Kirill Terekhov committed
350
351
352
353
354
355
356
			/// 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);}
Kirill Terekhov's avatar
Kirill Terekhov committed
357
358
359
360
361
362

			void                    Print() 
			{
				for(iterator it = Begin(); it != End(); ++it) std::cout << "(" << it->first << "," << it->second << ") ";
				std::cout << std::endl;
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
363
		};
Kirill Terekhov's avatar
Kirill Terekhov committed
364
365
366


		
Kirill Terekhov's avatar
Kirill Terekhov committed
367
		
368
369
370
		/// 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
371
372
373
374
375
376
377
378
379
380
381
382
		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:
383
384
385
386
387
			/// 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
388
389
390
391
392
			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();
			
393
			/// Return reference to i-th Row of the matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
394
			Row & operator [](INMOST_DATA_ENUM_TYPE i) {return data[i];}
395
			/// Return reference to i-th Row of the matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
396
			const Row & operator [](INMOST_DATA_ENUM_TYPE i) const {return data[i];}
397
			/// Return the total number of rows in the matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
398
399
400
401
402
403
			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();}
404
			/// Set the start and the end row numbers of the distributed matrix interval.
Kirill Terekhov's avatar
Kirill Terekhov committed
405
			void                 SetInterval(INMOST_DATA_ENUM_TYPE   start, INMOST_DATA_ENUM_TYPE   end)       {data.set_interval_beg(start); data.set_interval_end(end);}
406
			/// Get the start and the end row numbers of the distributed matrix interval.
Kirill Terekhov's avatar
Kirill Terekhov committed
407
408
			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);}
409
			/// Get the first row index of the distributed matrix interval.
Kirill Terekhov's avatar
Kirill Terekhov committed
410
			INMOST_DATA_ENUM_TYPE  GetFirstIndex() const {return data.get_interval_beg();}
411
			/// Get the communicator which the matrix is associated with.
Kirill Terekhov's avatar
Kirill Terekhov committed
412
413
414
415
416
417
418
419
420
421
422
423
424
425
			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;
			}
			
426
427
428
			/// Matrix-vector product of the form: y = alpha*A*x + beta * y.
			/// @param y Input/output vector.
			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
429
			void MatVecTranspose(INMOST_DATA_REAL_TYPE alpha, Solver::Vector & x, INMOST_DATA_REAL_TYPE beta, Solver::Vector & y) const; //y = alpha*At*x + beta * y
Kirill Terekhov's avatar
Kirill Terekhov committed
430

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

434
435
436
			/// 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
437
			void				 Load(std::string file, INMOST_DATA_ENUM_TYPE beg = ENUMUNDEF, INMOST_DATA_ENUM_TYPE end = ENUMUNDEF);
438
439
			/// 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
440
441
442
			void                 Save(std::string file);
			bool &               isParallel() { return is_parallel; }
			
443
			/// Get the matrix name specified in the main constructor.
Kirill Terekhov's avatar
Kirill Terekhov committed
444
445
446
			std::string          GetName() {return name;}
			//~ friend class Solver;
		};
Kirill Terekhov's avatar
Kirill Terekhov committed
447
448
449
450
451
452
453
454

		/// This class may be used to sum multiple sparse rows.
		/// \warning
		/// In parallel column indices of the matrix may span wider then 
		/// local row indices, to prevent any problem you are currently
		/// advised to set total size of the matrix as interval of the
		/// RowMerger. In future this may change, see todo 2 below.
		/// \todo
Kirill Terekhov's avatar
Kirill Terekhov committed
455
		/// 1. (testing!) Add iterators over entries.
Kirill Terekhov's avatar
Kirill Terekhov committed
456
457
458
459
460
461
462
463
		/// 2. Implement multiple intervals for distributed computation,
		///    then in parallel the user may specify additional range of indexes
		///    for elements that lay on the borders between each pair of processors.
		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.
Kirill Terekhov's avatar
Kirill Terekhov committed
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
      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;}
        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;
      };
Kirill Terekhov's avatar
Kirill Terekhov committed
488
489
490
491
492
		private:
			bool Sorted; ///< Contents of linked list should be sorted.
			INMOST_DATA_ENUM_TYPE Nonzeros; ///< Number of nonzero in linked list.
			interval< INMOST_DATA_ENUM_TYPE, Row::entry > LinkedList; ///< Storage for linked list.
		public:
493
			/// Default constructor without size specified.
Kirill Terekhov's avatar
Kirill Terekhov committed
494
495
496
497
			RowMerger();
			/// Constructor with size specified.
			/// @param interval_begin First index in linked list.
			/// @param interval_end Last index in linked list.
498
			/// @param Sorted Result should be sorted or not.
Kirill Terekhov's avatar
Kirill Terekhov committed
499
			RowMerger(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool Sorted = true);
500
			/// Constructor that gets sizes from the matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
501
502
503
504
505
506
507
508
509
510
511
			/// @param A Matrix to get sizes from.
			/// @param Sorted Result should be sorted.
			RowMerger(Matrix & A, bool Sorted = true);
			/// Destructor.
			~RowMerger();
			/// Resize linked list for new interval.
			/// \warning
			/// All contents of linked list will be lost after resize.
			/// This behavior may be changed in future.
			/// @param interval_begin First index in linked list.
			/// @param interval_end Last index in linked list.
512
			/// @param Sorted Result should be sorted or not.
Kirill Terekhov's avatar
Kirill Terekhov committed
513
514
515
516
517
518
			void Resize(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool Sorted = true);
			/// Resize linked list for new matrix.
			/// \warning
			/// All contents of linked list will be lost after resize.
			/// This behavior may be changed in future.
			/// @param A Matrix to get sizes from.
519
			/// @param Sorted Result should be sorted or not.
Kirill Terekhov's avatar
Kirill Terekhov committed
520
521
522
523
			void Resize(Matrix & A, bool Sorted = true);
			/// Clear linked list.
			void Clear();
			/// Add a row with a coefficient into empty linked list.
Kirill Terekhov's avatar
Kirill Terekhov committed
524
			/// This routine should be a bit faster then Solver::RowMerger::AddRow
Kirill Terekhov's avatar
Kirill Terekhov committed
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
			/// 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 = false);
			/// Add a row with a coefficient into non-empty linked list.
			/// Use Solver::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 = false);
			/// Multiply all entries of linked list by a coefficient.
			/// @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;}
Kirill Terekhov's avatar
Kirill Terekhov committed
550
551
552
553
554
555
556
557
558
559
560
      /// 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;

      iterator Begin() {return iterator(&LinkedList);}
      iterator End() {iterator ret(&LinkedList); ret.pos = EOL; return ret;}
Kirill Terekhov's avatar
Kirill Terekhov committed
561
		};
Kirill Terekhov's avatar
Kirill Terekhov committed
562
563
		
	private:
Kirill Terekhov's avatar
Kirill Terekhov committed
564
		static bool is_initialized, is_finalized;
Kirill Terekhov's avatar
Kirill Terekhov committed
565
566
567
568
569
570
		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;
Kirill Terekhov's avatar
Kirill Terekhov committed
571
572
573
574
575
576

		INMOST_DATA_ENUM_TYPE additive_schwartz_overlap;

		INMOST_DATA_ENUM_TYPE maximum_iterations;
		INMOST_DATA_REAL_TYPE absolute_tolerance;
		INMOST_DATA_REAL_TYPE relative_tolerance;
577
		INMOST_DATA_REAL_TYPE divergence_tolerance;
Kirill Terekhov's avatar
Kirill Terekhov committed
578
579
580
581
582

		INMOST_DATA_REAL_TYPE preconditioner_drop_tolerance;
		INMOST_DATA_REAL_TYPE preconditioner_reuse_tolerance;
		INMOST_DATA_REAL_TYPE preconditioner_ddpq_tolerance;
		INMOST_DATA_ENUM_TYPE preconditioner_reorder_nonzero;
Kirill Terekhov's avatar
Fix    
Kirill Terekhov committed
583
		INMOST_DATA_REAL_TYPE preconditioner_fill_level;
Kirill Terekhov's avatar
Kirill Terekhov committed
584
585
586
587
588
589
590
		INMOST_DATA_ENUM_TYPE preconditioner_rescale_iterations;
		INMOST_DATA_ENUM_TYPE preconditioner_condition_estimation;
		INMOST_DATA_ENUM_TYPE preconditioner_adapt_ddpq_tolerance;

		INMOST_DATA_ENUM_TYPE solver_gmres_substeps;

		std::string return_reason;
Kirill Terekhov's avatar
Kirill Terekhov committed
591
592
593
594
595
596
597
598
599
600
601
602
		
		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:
Kirill Terekhov's avatar
Kirill Terekhov committed
603
604
605
606
607
608
    /// Retrive approximate condition number produced by INNER_MPTILUC.
    /// The number is cond(L^-1).
    INMOST_DATA_REAL_TYPE GetConditionNumberL();
    /// Retrive approximate condition number produced by INNER_MPTILUC.
    /// The number is cond(U^-1).
    INMOST_DATA_REAL_TYPE GetConditionNumberU();
609
		/// Set the solver parameter of the integer type.
610
611
		/// You can find defaults for parameters in the top of the file inmost_solver.h.
		///
Kirill Terekhov's avatar
Kirill Terekhov committed
612
		/// Parameters:
613
614
		/// - "maximum_iterations" - total number of iterations
		/// - "schwartz_overlap"   - number of overlapping levels for additive schwartz method,
Kirill Terekhov's avatar
Kirill Terekhov committed
615
616
617
618
		///                          works for: 
		///                          INNER_ILU2, INNER_MLILUC
		///                          Trilinos_Aztec, Trilinos_Belos, Trilinos_ML, Trilinos_Ifpack
		///                          PETSc
619
		/// - "gmres_substeps"     - number of gmres steps performed after each bicgstab step,
Kirill Terekhov's avatar
Kirill Terekhov committed
620
621
		///                          works for:
		///                          INNER_ILU2, INNER_MLILUC
622
		/// - "reorder_nonzeros"   - place sparser rows at the beggining of matrix during reordering,
Kirill Terekhov's avatar
Kirill Terekhov committed
623
624
		///                          works for:
		///                          INNER_MLILUC
625
		/// - "rescale_iterations" - number of iterations for two-side matrix rescaling,
Kirill Terekhov's avatar
Kirill Terekhov committed
626
627
		///                          works for:
		///                          INNER_ILU2, INNER_MLILUC
628
629
		/// - "condition_estimation" - exploit condition estimation of inversed factors to adapt
		///                          drop and reuse tolerances,
Kirill Terekhov's avatar
Kirill Terekhov committed
630
631
		///                          works for:
		///                          INNER_MLILUC
632
		/// - "adapt_ddpq_tolerance" - adapt ddpq tolerance depending from the complexity 
633
		///                          of calculation of Schur complement,
Kirill Terekhov's avatar
Kirill Terekhov committed
634
635
		///                          works for:
		///                          INNER_MLILUC
Kirill Terekhov's avatar
Kirill Terekhov committed
636
		void SetParameterEnum(std::string name, INMOST_DATA_ENUM_TYPE value);
637
		/// Set the solver parameter of the real type.
638
639
		/// You can find defaults for parameters in the top of the file inmost_solver.h.
		///
Kirill Terekhov's avatar
Kirill Terekhov committed
640
		/// Parameters:
641
		/// - "absolute_tolerance" - iterative method will stop on i-th iteration
Kirill Terekhov's avatar
Kirill Terekhov committed
642
		///                          if ||A x(i)-b|| < absolute_tolerance
643
		/// - "relative_tolerance" - iterative method will stop on i-th iteration
Kirill Terekhov's avatar
Kirill Terekhov committed
644
		///                          if ||A x(i)-b||/||A x(0) - b||
645
		/// - "divergence_tolerance" - iterative method will fail if
Kirill Terekhov's avatar
Kirill Terekhov committed
646
		///                          ||A x(i) - b|| > divergence_tolerance
647
		/// - "drop_tolerance"     - tolerance for dropping values during incomplete factorization,
Kirill Terekhov's avatar
Kirill Terekhov committed
648
649
		///                          works for:
		///                          INNER_ILU2, INNER_MLILUC
Kirill Terekhov's avatar
Fix    
Kirill Terekhov committed
650
		///                          Trilinos_Aztec, Trilinos_Ifpack
Kirill Terekhov's avatar
Kirill Terekhov committed
651
		///                          PETSc
652
		/// - "reuse_tolerance"    - tolerance for reusing values during incomplete factorization,
Kirill Terekhov's avatar
Kirill Terekhov committed
653
		///                          these values are used only during calculation of L and U factors
654
655
656
		///                          and/or Schur complement and discarded once factorization is done,
		///                          value should be less then "drop_tolerance",
		///                          typical value is drop_tolerance^2,
Kirill Terekhov's avatar
Kirill Terekhov committed
657
658
		///                          works for:
		///                          INNER_ILU2, INNER_MLILUC
659
		/// - "ddpq_tolerance"     - by this tolerance most diagonnaly-dominant elements will be selected
Kirill Terekhov's avatar
Kirill Terekhov committed
660
661
662
663
664
665
666
		///                          to form the next level of factorization, the closer the tolerance
		///                          is to one the smaller will be the level. Actual rule is:
		///                          A(i,j)/(sum(A(i,:))+sum(A(:,j))-A(i,j)) > ddpq_tolerance *
		///                          A(imax,jmax)/(sum(A(imax,:))+sum(A(:,jmax))-A(imax,jmax))
		///                          where on imax, jmax maximum is reached.
		///                          works for:
		///                          INNER_MLILUC
667
		/// - "fill_level"         - level of fill for ILU-type preconditioners,
Kirill Terekhov's avatar
Fix    
Kirill Terekhov committed
668
669
670
		///                          works for:
		///                          INNER_ILU2 (if LFILL is defined in solver_ilu2.hpp)
		///                          Trilinos, Trilinos_Ifpack
Kirill Terekhov's avatar
Kirill Terekhov committed
671
		void SetParameterReal(std::string name, INMOST_DATA_REAL_TYPE value);
672
		/// Get the used defined name of the Solver.
Kirill Terekhov's avatar
Kirill Terekhov committed
673
674
		std::string          GetName() {return name;}
	
675
		/// Get the package Type.
Kirill Terekhov's avatar
Kirill Terekhov committed
676
677
		Type GetPackage() const {return _pack;}
		
678
		/// Set the matrix and construct the preconditioner.
679
		/// @param A Matrix A in linear problem Ax = b
680
681
682
		/// @param OldPreconditioner If this parameter is set to true,
		/// then the previous preconditioner will be used,
		/// otherwise the new preconditioner will be constructed. 
683
		///
684
		/// Preconditioner will be constructed on call to this function
685
686
687
688
		/// - for INNER_*, PETSc and ANI packages
		/// - for Trilinos preconditioner will be constructed each time Solver::Solve is called
		///
		/// Any changes to preconditioner parameters should happen before that point.
689
		/// If you increase gmres_substep after this point, inner methods most likely will fail
Kirill Terekhov's avatar
Kirill Terekhov committed
690
		void SetMatrix(Matrix & A, bool OldPreconditioner = false);
691
		/// Solver the linear system: A*x = b.
692
		/// Prior to this call you should call SetMatrix
693
		///
694
695
		/// @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.
696
		///
697
698
		/// It is assumed that the coefficient matrix A have been set
		/// and the preconditioner have been already constructed.
699
		///
700
		/// @see Solver::SetMatrix
Kirill Terekhov's avatar
Kirill Terekhov committed
701
		bool Solve(Vector & RHS, Vector & SOL);
702
703
704
705
706
707
		/// Return the number of iterations performed by the last solution.
		/// @see Solver::Solve
		INMOST_DATA_ENUM_TYPE Iterations();
		/// Return the final residual achieved by the last solution.
		/// @see Solver::Solve
		INMOST_DATA_REAL_TYPE Residual();
708
709
		/// Get the reason of convergence or divergence of the last solution.
		/// @see Solver::Solve
Kirill Terekhov's avatar
Kirill Terekhov committed
710
		std::string GetReason();
Kirill Terekhov's avatar
Kirill Terekhov committed
711

712
		/// Main constructor of the solver.
713
714
		/// Solver name provided here is used to extract options from database file
		/// for PETSc and Trilinos packages.
715
716
		/// @param pack The package Type to be used for solution.
		/// @param _name The user specified name of the current solver.
717
		/// @param comm Communicator for parallel data exchanges, MPI_COMM_WORLD by default.
718
719
720
721
		/// @see Solver::Initialize
		/// @see Solver::SetMatrix
		/// @see Solver::Solve
		/// @see Solver::Finalize
Kirill Terekhov's avatar
Kirill Terekhov committed
722
723
		Solver(Type pack, std::string _name = "", INMOST_MPI_Comm comm = INMOST_MPI_COMM_WORLD);
		~Solver();
724
725
		/// Initialize the stage of parallel solution.
		/// If MPI is not initialized yet, then it will be initialized.
726
		///
727
728
729
730
731
		/// database file is used to pass parameters to PETSc and Trilinos packages.
		/// if database file for is provided any changes through SetParameterEnum,
		/// SetParameterReal would not be effective for PETSc and Trilinos packages.
		/// Currently this database file provides directions for package-specific
		/// files. In future it is supposed to set up parameters for internal solvers.
732
733
734
		/// @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.
735
		///
736
737
738
		/// The shortest call to this function with the default solver parameters is the following: Initialize(NULL,NULL,"");
		/// @see Solver::Finalize
		/// @see Solver::isInitialized
739
		///
740
		/// Example of contents of the database file:
741
742
743
744
745
746
		///
		/// 	PETSc: petsc_options.txt
		/// 	Trilinos_Ifpack: trilinos_ifpack_options.xml
		/// 	Trilinos_ML: trilinos_ml_options.xml
		/// 	Trilinos_Aztec: trilinos_aztec_options.xml
		/// 	Trilinos_Belos: trilinos_belos_options.xml
Kirill Terekhov's avatar
Kirill Terekhov committed
747
		static void Initialize(int * argc, char *** argv, const char * database = "");
748
749
750
751
752
		/// 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
753
		static void Finalize();
Kirill Terekhov's avatar
Kirill Terekhov committed
754
755
756
		static bool isInitialized() {return is_initialized;}
		static bool isFinalized() {return is_finalized;}

757
		/// Clear all internal data of the current solver including matrix, preconditioner etc.
Kirill Terekhov's avatar
Kirill Terekhov committed
758
		void Clear();
Kirill Terekhov's avatar
Kirill Terekhov committed
759
760
761
762
763
764
	};
}

#endif // USE_SOLVER

#endif // INMOST_SOLVER_INCLUDED