inmost_solver.h 13.5 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1
2
3
#ifndef INMOST_SOLVER_INCLUDED
#define INMOST_SOLVER_INCLUDED

Dmitry Bagaev's avatar
Dmitry Bagaev committed
4
#include "inmost_solver_interface.h"
Kirill Terekhov's avatar
Kirill Terekhov committed
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#include "inmost_common.h"
#include "inmost_sparse.h"

#define DEFAULT_ADDITIVE_SCHWARTZ_OVERLAP             1
#define DEFAULT_ABSOLUTE_TOLERANCE                    1.0e-5
#define DEFAULT_RELATIVE_TOLERANCE                    1.0e-12
#define DEFAULT_DIVERGENCE_TOLERANCE                  1.0e+100
#define DEFAULT_MAXIMUM_ITERATIONS                    2500
#define DEFAULT_SOLVER_GMRES_SUBSTEPS                 2
#define DEFAULT_PRECONDITIONER_DROP_TOLERANCE         0.005
#define DEFAULT_PRECONDITIONER_REUSE_TOLERANCE        0.00005
#define DEFAULT_PRECONDITIONER_FILL_LEVEL             3
#define DEFAULT_PRECONDITIONER_DDPQ_TOLERANCE         0.75
#define DEFAULT_PRECONDITIONER_REORDER_NONZEROS       1
#define DEFAULT_PRECONDITIONER_RESCALE_ITERS          6
#define DEFAULT_PRECONDITIONER_CONDITION_ESTIMATION   1
#define DEFAULT_PRECONDITIONER_ADAPT_DDPQ_TOLERANCE   1

#if defined(USE_SOLVER)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
24
25
namespace INMOST
{
Dmitry Bagaev's avatar
Dmitry Bagaev committed
26
27
28
29

#define GUARD_MPI(x) {ierr = x; if( ierr != MPI_SUCCESS ) {char str[4096]; int len; MPI_Error_string(ierr,str,&len); std::cout << #x << " not successfull: " << str << std::endl; MPI_Abort(comm,-1000);}}
#define HASH_TABLE_SIZE 2048

Dmitry Bagaev's avatar
Dmitry Bagaev committed
30
31
    class Solver
    {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
32
33
34
35
36
37
38
39
40
41
42
43
44
    private:
        static int *argc;
        static char ***argv;
        static const char *database;
        static bool is_initialized;
        static bool is_finalized;

        //Actual solver using for solving system
        SolverInterface *solver;
        std::string prefix;
    public:

        /// Base class for low level operations with objects of Solver class.
Dmitry Bagaev's avatar
Dmitry Bagaev committed
45
46
        class OrderInfo
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
        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;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
63
            int rank, size;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
64
65
        public:
            void Clear();
Dmitry Bagaev's avatar
Dmitry Bagaev committed
66

Dmitry Bagaev's avatar
Dmitry Bagaev committed
67
            /// Return true if Matrix data have already been specified.
Dmitry Bagaev's avatar
Dmitry Bagaev committed
68
69
            bool &HaveMatrix() { return have_matrix; }

Dmitry Bagaev's avatar
Dmitry Bagaev committed
70
            OrderInfo();
Dmitry Bagaev's avatar
Dmitry Bagaev committed
71
72
73
74
75

            OrderInfo(const OrderInfo &other);

            OrderInfo &operator=(OrderInfo const &other);

Dmitry Bagaev's avatar
Dmitry Bagaev committed
76
            ~OrderInfo();
Dmitry Bagaev's avatar
Dmitry Bagaev committed
77

Dmitry Bagaev's avatar
Dmitry Bagaev committed
78
79
80
81
82
            /// 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.
Dmitry Bagaev's avatar
Dmitry Bagaev committed
83
84
            void PrepareMatrix(Sparse::Matrix &m, INMOST_DATA_ENUM_TYPE overlap);

Dmitry Bagaev's avatar
Dmitry Bagaev committed
85
            /// Restore initial nonparallel state of the Matrix with no overlap.
Dmitry Bagaev's avatar
Dmitry Bagaev committed
86
87
            void RestoreMatrix(Sparse::Matrix &m);

Dmitry Bagaev's avatar
Dmitry Bagaev committed
88
            /// Prepare parallel state of the Vector.
Dmitry Bagaev's avatar
Dmitry Bagaev committed
89
90
            void PrepareVector(Sparse::Vector &v) const;

Dmitry Bagaev's avatar
Dmitry Bagaev committed
91
            /// Restore initial nonparallel state of the Vector.
Dmitry Bagaev's avatar
Dmitry Bagaev committed
92
            void RestoreVector(Sparse::Vector &v) const;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
93
            /// Retrieve the processor number by binary search for the specified global index.
Dmitry Bagaev's avatar
Dmitry Bagaev committed
94
95
96
97
98
            INMOST_DATA_ENUM_TYPE GetProcessor(
                    INMOST_DATA_ENUM_TYPE gind) const; //retrieve 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;

Dmitry Bagaev's avatar
Dmitry Bagaev committed
99
            /// Get the local index region for the specified process.
Dmitry Bagaev's avatar
Dmitry Bagaev committed
100
101
102
            void
            GetLocalRegion(INMOST_DATA_ENUM_TYPE proc, INMOST_DATA_ENUM_TYPE &mbeg, INMOST_DATA_ENUM_TYPE &mend) const;

Dmitry Bagaev's avatar
Dmitry Bagaev committed
103
            /// Get the local index region for the current process.
Dmitry Bagaev's avatar
Dmitry Bagaev committed
104
105
            void GetVectorRegion(INMOST_DATA_ENUM_TYPE &mbeg, INMOST_DATA_ENUM_TYPE &mend) const
            {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
106
107
108
                mbeg = local_vector_begin;
                mend = local_vector_end;
            }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
109
            /// Get the rank of the current communicator, i.e. the current process index.
Dmitry Bagaev's avatar
Dmitry Bagaev committed
110
            INMOST_DATA_ENUM_TYPE GetRank() const { return rank; }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
111
            /// Get the size of the current communicator, i.e. the total number of processes used.
Dmitry Bagaev's avatar
Dmitry Bagaev committed
112
113
            INMOST_DATA_ENUM_TYPE GetSize() const { return size; }

Dmitry Bagaev's avatar
Dmitry Bagaev committed
114
            /// Update the shared data in parallel vector.
Dmitry Bagaev's avatar
Dmitry Bagaev committed
115
            void Update(Sparse::Vector &x); // update parallel vector
Dmitry Bagaev's avatar
Dmitry Bagaev committed
116
            /// Sum shared values in parallel vector.
Dmitry Bagaev's avatar
Dmitry Bagaev committed
117
            void Accumulate(Sparse::Vector &x); // sum shared values in parallel vector
Dmitry Bagaev's avatar
Dmitry Bagaev committed
118
            /// Get the sum of num elements of real array on all processes.
Dmitry Bagaev's avatar
Dmitry Bagaev committed
119
            void Integrate(INMOST_DATA_REAL_TYPE *inout, INMOST_DATA_ENUM_TYPE num) const;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
120
            /// Get the communicator which the solver is associated with.
Dmitry Bagaev's avatar
Dmitry Bagaev committed
121
            INMOST_MPI_Comm GetComm() const { return comm; }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
122
            // Access to arrays below allows to organize manual exchange
Dmitry Bagaev's avatar
Dmitry Bagaev committed
123
124
            INMOST_MPI_Request *GetSendRequests()
            {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
125
126
127
128
                assert(!send_requests.empty());
                return &send_requests[0];
            }

Dmitry Bagaev's avatar
Dmitry Bagaev committed
129
130
            INMOST_MPI_Request *GetRecvRequests()
            {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
131
132
133
134
135
136
137
138
                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()); }

Dmitry Bagaev's avatar
Dmitry Bagaev committed
139
140
            INMOST_DATA_ENUM_TYPE *GetSendExchangeArray()
            {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
141
142
143
144
145
146
                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()); }

Dmitry Bagaev's avatar
Dmitry Bagaev committed
147
148
            INMOST_DATA_ENUM_TYPE *GetRecvExchangeArray()
            {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
149
150
151
152
153
                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()); }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
            //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);}

            // 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;
        };

        // Main constructor of the solver.
        /// @param solverName The solver name to be used for solution.
        /// @param prefix 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
        Solver(std::string solverName, std::string prefix = "", INMOST_MPI_Comm _comm = INMOST_MPI_COMM_WORLD);

Dmitry Bagaev's avatar
Dmitry Bagaev committed
173
174
175
        Solver(const Solver &other);

        Solver &operator=(const Solver &other);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
176
177
178
179

        /// Return the solver name 
        /// @see Sparse::Solve
        std::string SolverName() const;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
180

Dmitry Bagaev's avatar
Dmitry Bagaev committed
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
        /// Return the solver user specified name of the current solver
        /// @see Sparse::Solve
        std::string SolverPrefix() const;

        /// Initialize the stage of parallel solution.
        /// If MPI is not initialized yet, then it will be initialized.
        ///
        /// 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.
        /// @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
        ///
        /// Example of contents of the database file:
        ///
        ///     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
        static void Initialize(int *argc, char ***argv, const char *database);

        /// 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
        static void Finalize();

        static bool isInitialized();
Dmitry Bagaev's avatar
Dmitry Bagaev committed
218

Dmitry Bagaev's avatar
Dmitry Bagaev committed
219
        static bool isFinalized();
Dmitry Bagaev's avatar
Dmitry Bagaev committed
220

Dmitry Bagaev's avatar
Dmitry Bagaev committed
221
222
223
224
225
226
227
228
229
230
231
232
233
234
        /// Set the matrix and construct the preconditioner.
        /// @param A Matrix A in linear problem Ax = b
        /// @param ModifiedPattern Indicates whether the structure of the matrix have 
        /// changed since last call to Solver::SetMatrix.
        /// @param OldPreconditioner If this parameter is set to true,
        /// then the previous preconditioner will be used,
        /// otherwise the new preconditioner will be constructed. 
        ///
        /// Preconditioner will be constructed on call to this function
        /// - for INNER_*, PETSc and ANI packages
        /// - for Trilinos preconditioner will be constructed each time Sparse::Solve is called
        ///
        /// Any changes to preconditioner parameters should happen before that point.
        /// If you increase gmres_substep after this point, inner methods most likely will fail
Dmitry Bagaev's avatar
Dmitry Bagaev committed
235
236
        void SetMatrix(Sparse::Matrix &A, bool ModifiedPattern = true, bool OldPreconditioner = false);

Dmitry Bagaev's avatar
Dmitry Bagaev committed
237
238
239
240
241
242
243
244
245
246
        /// Solver the linear system: A*x = b.
        /// Prior to this call you should call SetMatrix
        ///
        /// @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 Sparse::SetMatrix
Dmitry Bagaev's avatar
Dmitry Bagaev committed
247
        bool Solve(INMOST::Sparse::Vector &RHS, INMOST::Sparse::Vector &SOL);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
248
249
250
251

        /// Clear all internal data of the current solver including matrix, preconditioner etc.
        bool Clear();

Dmitry Bagaev's avatar
Dmitry Bagaev committed
252
        SolverParameters &GetParameters() const;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
253
254
255
256

        /// Return the number of iterations performed by the last solution.
        /// @see Sparse::Solve
        const INMOST_DATA_ENUM_TYPE Iterations() const;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
257

Dmitry Bagaev's avatar
Dmitry Bagaev committed
258
259
260
        /// Return the final residual achieved by the last solution.
        /// @see Sparse::Solve
        const INMOST_DATA_REAL_TYPE Residual() const;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
261

Dmitry Bagaev's avatar
Dmitry Bagaev committed
262
263
264
265
        /// Get the reason of convergence or divergence of the last solution.
        /// @see Sparse::Solve
        const std::string ReturnReason() const;

Dmitry Bagaev's avatar
Dmitry Bagaev committed
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
        /// Computes the smallest and the largest eigenvalue with the power method.
        /// Requires SetMatrix to be called to compute the preconditioner.
        /// Currently works for internal methods only, since it uses internal matrix-vector multiplication.
        /// Largest eigenvalue: vprev = 0; v = rand(); while( |v|-|vprev| > tol ) {vprev = v; v = A*v; v /= |v|;}
        ///                     lambda_max = |v|;
        /// Smallest eigenvalue: vprev = 0; v = rand(); while( |v|-|vprev| > tol ){vprev = v; solve(A*v = v); v /= |v|;}
        ///                     lambda_min = 1.0/|v|;
        /// See answer by Blair Perot in:
        /// https://www.researchgate.net/post/What_is_the_best_way_to_estimate_the_condition_number_of_a_sparse_matrix.
        /// @param tol Tolerance used for power series.
        /// @param maxits Maximum number of iterations allowed.
        /// @return Condition number or 1.0e100 if not converged.
        INMOST_DATA_REAL_TYPE Condest(INMOST_DATA_REAL_TYPE tol, INMOST_DATA_ENUM_TYPE maxits = 100);

        static bool isSolverAvailable(std::string name);

        static std::vector<std::string> getAvailableSolvers();

Dmitry Bagaev's avatar
Dmitry Bagaev committed
284
        ~Solver();
Dmitry Bagaev's avatar
Dmitry Bagaev committed
285
286

    private:
Dmitry Bagaev's avatar
Dmitry Bagaev committed
287
288
        static std::string parseDatabase(std::string solverName);
    };
Dmitry Bagaev's avatar
Dmitry Bagaev committed
289
290

    typedef std::vector<std::string>::iterator solvers_names_iterator_t;
Kirill Terekhov's avatar
Kirill Terekhov committed
291
292
}

Dmitry Bagaev's avatar
Dmitry Bagaev committed
293
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
294

Dmitry Bagaev's avatar
Dmitry Bagaev committed
295
#endif //INMOST_SOLVER_INCLUDED