SolverFCBIILU2.cpp 8.34 KB
Newer Older
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1
#include <Source/Misc/utils.h>
2
#include "SolverFCBIILU2.h"
Dmitry Bagaev's avatar
Dmitry Bagaev committed
3
#include "solver_fcbiilu2.h"
4
5
6

namespace INMOST {

Dmitry Bagaev's avatar
Dmitry Bagaev committed
7
    SolverFCBIILU2::SolverFCBIILU2() {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
8
9
        solver_data = NULL;
        matrix_data = NULL;
10
11
    }

Dmitry Bagaev's avatar
Dmitry Bagaev committed
12
    SolverInterface *SolverFCBIILU2::Copy(const SolverInterface *other) {
13
14
15
16
17
18
19
20
21
        if (other == NULL) {
            throw INMOST::SolverCopyNullException;
        }
        const SolverFCBIILU2 *fcother;
        try {
            fcother = dynamic_cast<const SolverFCBIILU2 *>(other);
        } catch (...) {
            throw INMOST::SolverCopyException;
        }
22
23
24
25
26
        SolverCopyDataFcbiilu2(&solver_data, fcother->solver_data, communicator);
        if (fcother->matrix_data != NULL) {
            MatrixCopyDataFcbiilu2(&matrix_data, fcother->matrix_data);
            SolverSetMatrixFcbiilu2(solver_data, matrix_data, false, false);
        }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
27
        return this;
28
29
30
    }

    void SolverFCBIILU2::Assign(const SolverInterface *other) {
31
32
33
34
35
36
37
38
39
        if (other == NULL) {
            throw INMOST::SolverAssignNullException;
        }
        const SolverFCBIILU2 *fcother;
        try {
            fcother = dynamic_cast<const SolverFCBIILU2 *>(other);
        } catch (...) {
            throw INMOST::SolverAssignException;
        }
40
41
42
43
44
45
46
47
48
49
50
        SolverAssignDataFcbiilu2(solver_data, fcother->solver_data);
        if (fcother->matrix_data != NULL) {
            if (matrix_data != NULL) {
                MatrixAssignDataFcbiilu2(matrix_data, fcother->matrix_data);
            } else {
                MatrixCopyDataFcbiilu2(&matrix_data, fcother->matrix_data);
            }
            SolverSetMatrixFcbiilu2(solver_data, matrix_data, false, false);
        }
    }

Dmitry Bagaev's avatar
Dmitry Bagaev committed
51
52
    void SolverFCBIILU2::Setup(int *argc, char ***argv, SolverParameters &p) {
        SolverInitDataFcbiilu2(&solver_data, communicator, p.solverPrefix.c_str());
Dmitry Bagaev's avatar
Dmitry Bagaev committed
53
54
55
56
57
        solver_data->kovl = 0;    // number of overlap layers: kovl=0,1,2,...
        solver_data->tau = 3e-3;  // the ILU2 precision (for the submatrix factorization); tau=3e-3
        solver_data->eps = 1e-5;  // the residual precision: ||r|| < eps * ||b||; eps=1e-6
        solver_data->nit = 999;   // number of iterations permitted; nit=999
        solver_data->msglev = 2;  // messages level; msglev=0 for silent; msglev=1 to output solution statistics
Dmitry Bagaev's avatar
Dmitry Bagaev committed
58
59
60
61
62
63
64
        if (p.internalFile != "") {
            SolverInitializeFcbiilu2(solver_data, argc, argv, p.internalFile.c_str());
        } else {
            for (parameters_iterator_t parameter = p.parameters.begin(); parameter < p.parameters.end(); parameter++) {
                this->SetParameter((*parameter).first, (*parameter).second);
            }
        }
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
    }

    void SolverFCBIILU2::SetMatrix(Sparse::Matrix &A, bool ModifiedPattern, bool OldPreconditioner) {
        bool modified_pattern = ModifiedPattern;
        //~ if( A.comm != comm ) throw DifferentCommunicatorInSolver;
        if (matrix_data == NULL) {
            MatrixInitDataFcbiilu2(&matrix_data, A.GetCommunicator(), A.GetName().c_str());
            modified_pattern = true;
        }
        if (modified_pattern) {
            global_size = local_size = A.Size();

            MatrixDestroyDataFcbiilu2(&matrix_data);
            MatrixInitDataFcbiilu2(&matrix_data, A.GetCommunicator(), A.GetName().c_str());
            INMOST_DATA_ENUM_TYPE nnz = 0, k = 0, q = 1, shift = 1;
            int nproc, myid;
#if defined(USE_MPI)
            MPI_Comm_size(A.GetCommunicator(), &nproc);
            MPI_Comm_rank(A.GetCommunicator(), &myid);
#else
            nproc = 1;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
86
            myid = 0;
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
#endif
            int *ibl = (int *) malloc(sizeof(int) * (nproc + 1));
            ibl[0] = 0;
            int n = A.Size();
#if defined(USE_MPI)
            INMOST_DATA_ENUM_TYPE mbeg, mend;
            A.GetInterval(mbeg, mend);
            n = mend - mbeg;
            int block_end = mend;
            MPI_Allgather(&block_end, 1, MPI_INT, &ibl[1], 1, MPI_INT, A.GetCommunicator());
#else
            ibl[1] = n;
#endif
            int *ia = (int *) malloc(sizeof(int) * (n + 1));
            for (Sparse::Matrix::iterator it = A.Begin(); it != A.End(); it++) nnz += it->Size();
            int *ja = (int *) malloc(sizeof(int) * nnz);
            double *values = (double *) malloc(sizeof(double) * nnz);
            ia[0] = shift;
            for (Sparse::Matrix::iterator it = A.Begin(); it != A.End(); it++) {
                for (Sparse::Row::iterator jt = it->Begin(); jt != it->End(); jt++) {
                    ja[k] = jt->first + 1;
                    values[k] = jt->second;
                    //std::cout<<"# q="<<q<<" k="<<k<<" ja="<<ja[k]<<" a="<<values[k]<<std::endl;//db!
                    k++;
                }
                shift += it->Size();
                ia[q++] = shift;
            }
            MatrixFillFcbiilu2(matrix_data, n, nproc, ibl, ia, ja, values);
        } else {
            INMOST_DATA_ENUM_TYPE nnz = 0, k = 0;
            for (Sparse::Matrix::iterator it = A.Begin(); it != A.End(); it++) nnz += it->Size();
            double *values = (double *) malloc(sizeof(double) * nnz);
            k = 0;
            for (Sparse::Matrix::iterator it = A.Begin(); it != A.End(); it++)
                for (Sparse::Row::iterator jt = it->Begin(); jt != it->End(); jt++)
                    values[k++] = jt->second;
            MatrixFillValuesFcbiilu2(matrix_data, values);
        }
        MatrixFinalizeFcbiilu2(matrix_data);
        SolverSetMatrixFcbiilu2(solver_data, matrix_data, modified_pattern, OldPreconditioner);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
128
        time_prec = solver_data->dstat[7];
129
130
131
132
133
134
    }

    bool SolverFCBIILU2::Solve(INMOST::Sparse::Vector &RHS, INMOST::Sparse::Vector &SOL) {
        INMOST_DATA_ENUM_TYPE vbeg, vend;
        RHS.GetInterval(vbeg, vend);

Dmitry Bagaev's avatar
Dmitry Bagaev committed
135
        vector *rhs_data = NULL;
136
137
138
        VectorInitDataFcbiilu2(&rhs_data, RHS.GetCommunicator(), RHS.GetName().c_str());
        VectorPreallocateFcbiilu2(rhs_data, local_size);

Dmitry Bagaev's avatar
Dmitry Bagaev committed
139
        vector *solution_data = NULL;
140
141
142
143
144
145
146
147
148
        VectorInitDataFcbiilu2(&solution_data, SOL.GetCommunicator(), SOL.GetName().c_str());
        VectorPreallocateFcbiilu2(solution_data, local_size);
        VectorFillFcbiilu2(rhs_data, &RHS[vbeg]);
        VectorFinalizeFcbiilu2(rhs_data);

        VectorFillFcbiilu2(solution_data, &SOL[vbeg]);
        VectorFinalizeFcbiilu2(solution_data);
        bool result = SolverSolveFcbiilu2(solver_data, rhs_data, solution_data);
        if (result) VectorLoadFcbiilu2(solution_data, &SOL[vbeg]);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
149
        iter_time = solver_data->dstat[9];
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
        return result;
    }

    bool SolverFCBIILU2::Clear() {
        if (matrix_data != NULL) {
            MatrixDestroyDataFcbiilu2(&matrix_data);
        }
        SolverDestroyDataFcbiilu2(&solver_data);
        return true;
    }

    bool SolverFCBIILU2::isMatrixSet() {
        return matrix_data != NULL;
    }

Dmitry Bagaev's avatar
Dmitry Bagaev committed
165
    std::string SolverFCBIILU2::GetParameter(std::string name) const {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
166
167
168
169
170
171
        if (name == "time_prec") return INMOST::to_string(time_prec);
        if (name == "time_iter") return INMOST::to_string(iter_time);
        if (name == "time_total") return INMOST::to_string(time_prec + iter_time);
        if (name == "prec_density") return INMOST::to_string(solver_data->dstat[0]);
        if (name == "pivot_mod") return INMOST::to_string(solver_data->istat[0]);
        std::cout << "Parameter " << name << " is unknown" << std::endl;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
172
173
174
175
        return "";
    }

    void SolverFCBIILU2::SetParameter(std::string name, std::string value) {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
176
177
178
179
180
181
182
        const char *val = value.c_str();
        if (name == "kovl") solver_data->kovl = atoi(val);
        else if (name == "tau") solver_data->tau = atof(val);
        else if (name == "eps") solver_data->eps = atof(val);
        else if (name == "nit") solver_data->nit = atoi(val);
        else if (name == "msglev") solver_data->msglev = atoi(val);
        else std::cout << "Parameter " << name << " is unknown" << std::endl;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
183
184
    }

185
    const INMOST_DATA_ENUM_TYPE SolverFCBIILU2::Iterations() const {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
186
        return static_cast<INMOST_DATA_ENUM_TYPE>(solver_data->ITER);
187
188
189
    }

    const INMOST_DATA_REAL_TYPE SolverFCBIILU2::Residual() const {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
190
        return solver_data->RESID;
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
    }

    const std::string SolverFCBIILU2::ReturnReason() const {
        return "Unspecified for FCBIILU2";
    }

    const std::string SolverFCBIILU2::SolverName() const {
        return "fcbiilu2";
    }

    void SolverFCBIILU2::Finalize() {
        SolverFinalizeFcbiilu2();
    }

    SolverFCBIILU2::~SolverFCBIILU2() {
        this->Clear();
    }

}