inmost_solver.h 29.8 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 45
			INNER_ILU2,     ///< inner Solver based on BiCGStab(L) solver with second order ILU factorization as preconditioner.
			INNER_MLILUC,   ///< 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.
Kirill Terekhov's avatar
Kirill Terekhov committed
46 47 48 49
			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
50
			PETSc,          ///< external Solver PETSc, @see http://www.mcs.anl.gov/petsc/.
Kirill Terekhov's avatar
Kirill Terekhov committed
51
			ANI             ///< external Solver from ANI3D based on ILU2 (sequential Fortran version).
Kirill Terekhov's avatar
Kirill Terekhov committed
52 53 54 55 56 57 58 59 60
		};

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

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

			/// Get the scalar product of the specified interval of the distributed vector.
			INMOST_DATA_REAL_TYPE ScalarProd(Vector const & left, Vector const & right, INMOST_DATA_ENUM_TYPE index_begin, INMOST_DATA_ENUM_TYPE index_end) const;
Kirill Terekhov's avatar
Kirill Terekhov committed
135 136
		};

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

Kirill Terekhov's avatar
Kirill Terekhov committed
183
			
Kirill Terekhov's avatar
Kirill Terekhov committed
184 185


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

197
			/// Clear all data of the current vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
198
			void Clear() {data.clear();}
Kirill Terekhov's avatar
Kirill Terekhov committed
199 200 201 202 203
			//~ friend class Solver;
		};



204
		/// Class to store the sparse matrix row.
Kirill Terekhov's avatar
Kirill Terekhov committed
205 206 207 208
		class Row 
		{
		public:
			
209
			/// Entry of the sparse matrix row.
Kirill Terekhov's avatar
Kirill Terekhov committed
210 211
			typedef struct entry_s 
			{
212 213
				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
214 215 216 217
				//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
218 219
				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
220 221 222 223 224 225 226
			__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
227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245
			
			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
246
			void Report() {data.report_addr();}
Kirill Terekhov's avatar
Kirill Terekhov committed
247 248 249 250 251 252 253 254 255 256 257 258
			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
259 260 261
				//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
262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287
#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; }
288
			/// 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
289 290 291 292 293 294
			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
295
				for(Entries::size_type it = 0; it < data.size(); ++it)
Kirill Terekhov's avatar
Kirill Terekhov committed
296 297 298 299 300 301 302 303 304
					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;
				
			}
305
			/// The operator [] used to access individual elements of the row.
Kirill Terekhov's avatar
Kirill Terekhov committed
306 307 308 309 310
			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
311
				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
312 313 314 315 316 317

				//you should not come here
				assert(false);
				return 1.0e20;
			}
			//void           Reserve(INMOST_DATA_ENUM_TYPE num) { data.reserve(num);}
318
			/// Clear all data of the current row.
Kirill Terekhov's avatar
Kirill Terekhov committed
319 320
			void                  Clear() { data.clear(); }
			void                  Swap(Solver::Row & other) { data.swap(other.data); bool tmp = marker; marker = other.marker; other.marker = tmp; }
321
			/// The size of the sparse row, i.e. the total number of nonzero elements.
Kirill Terekhov's avatar
Kirill Terekhov committed
322 323 324 325 326 327 328 329 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;}
			
			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(); }
336
			/// Return the scalar product of the current sparse row by a dense Vector.
Kirill Terekhov's avatar
Kirill Terekhov committed
337 338
			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
339
			/// Set the vector entries by zeroes.
Kirill Terekhov's avatar
Kirill Terekhov committed
340
			void                  Zero() {for(iterator it = Begin(); it != End(); ++it) it->second = 0;}
341 342 343
			/// 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
Fixes  
Kirill Terekhov committed
344
			void                  Push(INMOST_DATA_ENUM_TYPE ind, INMOST_DATA_REAL_TYPE val) {data.push_back(make_entry(ind,val));}
Kirill Terekhov's avatar
Kirill Terekhov committed
345 346
		};
		
347 348 349
		/// 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
350 351 352 353 354 355 356 357 358 359 360 361
		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:
362 363 364 365 366
			/// 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
367 368 369 370 371
			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();
			
372
			/// Return reference to i-th Row of the matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
373
			Row & operator [](INMOST_DATA_ENUM_TYPE i) {return data[i];}
374
			/// Return reference to i-th Row of the matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
375
			const Row & operator [](INMOST_DATA_ENUM_TYPE i) const {return data[i];}
376
			/// Return the total number of rows in the matrix.
Kirill Terekhov's avatar
Kirill Terekhov committed
377 378 379 380 381 382
			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();}
383
			/// Set the start and the end row numbers of the distributed matrix interval.
Kirill Terekhov's avatar
Kirill Terekhov committed
384
			void                 SetInterval(INMOST_DATA_ENUM_TYPE   start, INMOST_DATA_ENUM_TYPE   end)       {data.set_interval_beg(start); data.set_interval_end(end);}
385
			/// Get the start and the end row numbers of the distributed matrix interval.
Kirill Terekhov's avatar
Kirill Terekhov committed
386 387
			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);}
388
			/// Get the first row index of the distributed matrix interval.
Kirill Terekhov's avatar
Kirill Terekhov committed
389
			INMOST_DATA_ENUM_TYPE  GetFirstIndex() const {return data.get_interval_beg();}
390
			/// Get the communicator which the matrix is associated with.
Kirill Terekhov's avatar
Kirill Terekhov committed
391 392 393 394 395 396 397 398 399 400 401 402 403 404
			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;
			}
			
405 406 407 408
			/// 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
409

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

413 414 415
			/// 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
416
			void				 Load(std::string file, INMOST_DATA_ENUM_TYPE beg = ENUMUNDEF, INMOST_DATA_ENUM_TYPE end = ENUMUNDEF);
417 418
			/// 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
419 420 421
			void                 Save(std::string file);
			bool &               isParallel() { return is_parallel; }
			
422
			/// Get the matrix name specified in the main constructor.
Kirill Terekhov's avatar
Kirill Terekhov committed
423 424 425 426 427
			std::string          GetName() {return name;}
			//~ friend class Solver;
		};
		
	private:
Kirill Terekhov's avatar
Kirill Terekhov committed
428
		static bool is_initialized, is_finalized;
Kirill Terekhov's avatar
Kirill Terekhov committed
429 430 431 432 433 434
		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
435 436 437 438 439 440

		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;
441
		INMOST_DATA_REAL_TYPE divergence_tolerance;
Kirill Terekhov's avatar
Kirill Terekhov committed
442 443 444 445 446

		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
447
		INMOST_DATA_REAL_TYPE preconditioner_fill_level;
Kirill Terekhov's avatar
Kirill Terekhov committed
448 449 450 451 452 453 454
		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
455 456 457 458 459 460 461 462 463 464 465 466
		
		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:
467
		/// Set the solver parameter of the integer type.
468 469
		/// You can find defaults for parameters in the top of the file inmost_solver.h.
		///
Kirill Terekhov's avatar
Kirill Terekhov committed
470
		/// Parameters:
471 472
		/// - "maximum_iterations" - total number of iterations
		/// - "schwartz_overlap"   - number of overlapping levels for additive schwartz method,
Kirill Terekhov's avatar
Kirill Terekhov committed
473 474 475 476
		///                          works for: 
		///                          INNER_ILU2, INNER_MLILUC
		///                          Trilinos_Aztec, Trilinos_Belos, Trilinos_ML, Trilinos_Ifpack
		///                          PETSc
477
		/// - "gmres_substeps"     - number of gmres steps performed after each bicgstab step,
Kirill Terekhov's avatar
Kirill Terekhov committed
478 479
		///                          works for:
		///                          INNER_ILU2, INNER_MLILUC
480
		/// - "reorder_nonzeros"   - place sparser rows at the beggining of matrix during reordering,
Kirill Terekhov's avatar
Kirill Terekhov committed
481 482
		///                          works for:
		///                          INNER_MLILUC
483
		/// - "rescale_iterations" - number of iterations for two-side matrix rescaling,
Kirill Terekhov's avatar
Kirill Terekhov committed
484 485
		///                          works for:
		///                          INNER_ILU2, INNER_MLILUC
486 487
		/// - "condition_estimation" - exploit condition estimation of inversed factors to adapt
		///                          drop and reuse tolerances,
Kirill Terekhov's avatar
Kirill Terekhov committed
488 489
		///                          works for:
		///                          INNER_MLILUC
490 491
		/// - "adapt_ddpq_tolerance" - adapt ddpq tolerance depending from the complexity 
		//                           of calculation of Schur complement,
Kirill Terekhov's avatar
Kirill Terekhov committed
492 493
		///                          works for:
		///                          INNER_MLILUC
Kirill Terekhov's avatar
Kirill Terekhov committed
494
		void SetParameterEnum(std::string name, INMOST_DATA_ENUM_TYPE value);
495
		/// Set the solver parameter of the real type.
496 497
		/// You can find defaults for parameters in the top of the file inmost_solver.h.
		///
Kirill Terekhov's avatar
Kirill Terekhov committed
498
		/// Parameters:
499
		/// - "absolute_tolerance" - iterative method will stop on i-th iteration
Kirill Terekhov's avatar
Kirill Terekhov committed
500
		///                          if ||A x(i)-b|| < absolute_tolerance
501
		/// - "relative_tolerance" - iterative method will stop on i-th iteration
Kirill Terekhov's avatar
Kirill Terekhov committed
502
		///                          if ||A x(i)-b||/||A x(0) - b||
503
		/// - "divergence_tolerance" - iterative method will fail if
Kirill Terekhov's avatar
Kirill Terekhov committed
504
		///                          ||A x(i) - b|| > divergence_tolerance
505
		/// - "drop_tolerance"     - tolerance for dropping values during incomplete factorization,
Kirill Terekhov's avatar
Kirill Terekhov committed
506 507
		///                          works for:
		///                          INNER_ILU2, INNER_MLILUC
Kirill Terekhov's avatar
Fix  
Kirill Terekhov committed
508
		///                          Trilinos_Aztec, Trilinos_Ifpack
Kirill Terekhov's avatar
Kirill Terekhov committed
509
		///                          PETSc
510
		/// - "reuse_tolerance"    - tolerance for reusing values during incomplete factorization,
Kirill Terekhov's avatar
Kirill Terekhov committed
511
		///                          these values are used only during calculation of L and U factors
512 513 514
		///                          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
515 516
		///                          works for:
		///                          INNER_ILU2, INNER_MLILUC
517
		/// - "ddpq_tolerance"     - by this tolerance most diagonnaly-dominant elements will be selected
Kirill Terekhov's avatar
Kirill Terekhov committed
518 519 520 521 522 523 524
		///                          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
525
		/// - "fill_level"         - level of fill for ILU-type preconditioners,
Kirill Terekhov's avatar
Fix  
Kirill Terekhov committed
526 527 528
		///                          works for:
		///                          INNER_ILU2 (if LFILL is defined in solver_ilu2.hpp)
		///                          Trilinos, Trilinos_Ifpack
Kirill Terekhov's avatar
Kirill Terekhov committed
529
		void SetParameterReal(std::string name, INMOST_DATA_REAL_TYPE value);
530
		/// Get the used defined name of the Solver.
Kirill Terekhov's avatar
Kirill Terekhov committed
531 532
		std::string          GetName() {return name;}
	
533
		/// Get the package Type.
Kirill Terekhov's avatar
Kirill Terekhov committed
534 535
		Type GetPackage() const {return _pack;}
		
536 537
		/// Return the number of iterations performed by the last solution.
		/// @see Solver::Solve
Kirill Terekhov's avatar
Kirill Terekhov committed
538
		INMOST_DATA_ENUM_TYPE Iterations();
539 540
		/// Return the final residual achieved by the last solution.
		/// @see Solver::Solve
Kirill Terekhov's avatar
Kirill Terekhov committed
541
		INMOST_DATA_REAL_TYPE Residual();
542
		/// Set the matrix and construct the preconditioner.
543
		/// @param A Matrix A in linear problem Ax = b
544 545 546
		/// @param OldPreconditioner If this parameter is set to true,
		/// then the previous preconditioner will be used,
		/// otherwise the new preconditioner will be constructed. 
547
		///
548
		/// Preconditioner will be constructed on call to this function
549 550 551 552
		/// - 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.
553
		/// If you increase gmres_substep after this point, inner methods most likely will fail
Kirill Terekhov's avatar
Kirill Terekhov committed
554
		void SetMatrix(Matrix & A, bool OldPreconditioner = false);
555
		/// Solver the linear system: A*x = b.
556
		/// Prior to this call you should call SetMatrix
557
		///
558 559
		/// @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.
560
		///
561 562
		/// It is assumed that the coefficient matrix A have been set
		/// and the preconditioner have been already constructed.
563
		///
564
		/// @see Solver::SetMatrix
Kirill Terekhov's avatar
Kirill Terekhov committed
565
		bool Solve(Vector & RHS, Vector & SOL);
566 567
		/// Get the reason of convergence or divergence of the last solution.
		/// @see Solver::Solve
Kirill Terekhov's avatar
Kirill Terekhov committed
568
		std::string GetReason();
Kirill Terekhov's avatar
Kirill Terekhov committed
569

570
		/// Main constructor of the solver.
571 572
		/// Solver name provided here is used to extract options from database file
		/// for PETSc and Trilinos packages.
573 574
		/// @param pack The package Type to be used for solution.
		/// @param _name The user specified name of the current solver.
575
		/// @param comm Communicator for parallel data exchanges, MPI_COMM_WORLD by default.
576 577 578 579
		/// @see Solver::Initialize
		/// @see Solver::SetMatrix
		/// @see Solver::Solve
		/// @see Solver::Finalize
Kirill Terekhov's avatar
Kirill Terekhov committed
580 581
		Solver(Type pack, std::string _name = "", INMOST_MPI_Comm comm = INMOST_MPI_COMM_WORLD);
		~Solver();
582 583
		/// Initialize the stage of parallel solution.
		/// If MPI is not initialized yet, then it will be initialized.
584
		///
585 586 587 588 589
		/// 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.
590 591 592
		/// @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.
593
		///
594 595 596
		/// The shortest call to this function with the default solver parameters is the following: Initialize(NULL,NULL,"");
		/// @see Solver::Finalize
		/// @see Solver::isInitialized
597
		///
598
		/// Example of contents of the database file:
599 600 601 602 603 604
		///
		/// 	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
605
		static void Initialize(int * argc, char *** argv, const char * database = "");
606 607 608 609 610
		/// 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
611
		static void Finalize();
Kirill Terekhov's avatar
Kirill Terekhov committed
612 613 614
		static bool isInitialized() {return is_initialized;}
		static bool isFinalized() {return is_finalized;}

615
		/// Clear all internal data of the current solver including matrix, preconditioner etc.
Kirill Terekhov's avatar
Kirill Terekhov committed
616
		void Clear();
Kirill Terekhov's avatar
Kirill Terekhov committed
617 618 619 620 621 622
	};
}

#endif // USE_SOLVER

#endif // INMOST_SOLVER_INCLUDED