SolverK3BIILU2.cpp 7.12 KB
Newer Older
1
#include "SolverK3BIILU2.h"
Dmitry Bagaev's avatar
Dmitry Bagaev committed
2
#include "solver_k3biilu2.h"
3
4
5

namespace INMOST {

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

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

    void SolverK3BIILU2::Assign(const SolverInterface *other) {
30
31
32
33
34
35
36
37
38
        if (other == NULL) {
            throw INMOST::SolverAssignNullException;
        }
        const SolverK3BIILU2 *k3other;
        try {
            k3other = dynamic_cast<const SolverK3BIILU2 *>(other);
        } catch (...) {
            throw INMOST::SolverAssignException;
        }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
39
40
41
42
43
44
45
46
47
48
49
        SolverAssignDataK3biilu2(solver_data, k3other->solver_data);
        if (k3other->matrix_data != NULL) {
            if (matrix_data != NULL) {
                MatrixAssignDataK3biilu2(matrix_data, k3other->matrix_data);
            } else {
                MatrixCopyDataK3biilu2(&matrix_data, k3other->matrix_data);
            }
            SolverSetMatrixK3biilu2(solver_data, matrix_data, false, false);
        }
    }

Dmitry Bagaev's avatar
Dmitry Bagaev committed
50
51
    void SolverK3BIILU2::Setup(int *argc, char ***argv, SolverParameters &p) {
        SolverInitDataK3biilu2(&solver_data, communicator, p.solverPrefix.c_str());
Dmitry Bagaev's avatar
Dmitry Bagaev committed
52
        SolverInitializeK3biilu2(solver_data, argc, argv, p.internalFile.c_str());
Dmitry Bagaev's avatar
Dmitry Bagaev committed
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
    }

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

            MatrixDestroyDataK3biilu2(&matrix_data);
            MatrixInitDataK3biilu2(&matrix_data, A.GetCommunicator(), A.GetName().c_str());
            INMOST_DATA_ENUM_TYPE nnz = 0, k = 0, q = 1, shift = 0;
            int nproc, myid;
69
#if defined(USE_MPI)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
70
71
            MPI_Comm_size(A.GetCommunicator(), &nproc);
            MPI_Comm_rank(A.GetCommunicator(), &myid);
72
#else
Dmitry Bagaev's avatar
Dmitry Bagaev committed
73
74
            nproc = 1;
            myid = 0;
75
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
76
77
78
            int *ibl = (int *) malloc(sizeof(int) * (nproc + 1));
            ibl[0] = 0;
            int n = A.Size();
79
#if defined(USE_MPI)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
80
81
82
83
84
            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());
85
#else
Dmitry Bagaev's avatar
Dmitry Bagaev committed
86
            ibl[1] = n;
87
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
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
            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 + 0;
                    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;
            }
            MatrixFillK3biilu2(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;
                }
            }
            MatrixFillValuesK3biilu2(matrix_data, values);
        }
        MatrixFinalizeK3biilu2(matrix_data);
        SolverSetMatrixK3biilu2(solver_data, matrix_data, modified_pattern, OldPreconditioner);
    }

    bool SolverK3BIILU2::Solve(INMOST::Sparse::Vector &RHS, INMOST::Sparse::Vector &SOL) {
        INMOST_DATA_ENUM_TYPE vbeg, vend;
126
        RHS.GetInterval(vbeg, vend);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
127
128
        vector_k3biilu2 *rhs_data = NULL;
        vector_k3biilu2 *solution_data = NULL;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165

        VectorInitDataK3biilu2(&rhs_data, RHS.GetCommunicator(), RHS.GetName().c_str());
        VectorPreallocateK3biilu2(rhs_data, local_size);

        VectorInitDataK3biilu2(&solution_data, SOL.GetCommunicator(), SOL.GetName().c_str());
        VectorPreallocateK3biilu2(solution_data, local_size);

        VectorFillK3biilu2(rhs_data, &RHS[vbeg]);
        VectorFinalizeK3biilu2(rhs_data);

        VectorFillK3biilu2(solution_data, &SOL[vbeg]);
        VectorFinalizeK3biilu2(solution_data);

        bool result = SolverSolveK3biilu2(solver_data, rhs_data, solution_data);
        if (result) VectorLoadK3biilu2(solution_data, &SOL[vbeg]);
        return result;
    }

    bool SolverK3BIILU2::Clear() {
        if (matrix_data != NULL) {
            MatrixDestroyDataK3biilu2(&matrix_data);
        }
        SolverDestroyDataK3biilu2(&solver_data);
        return true;
    }

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

    std::string SolverK3BIILU2::GetParameter(std::string name) const {
        std::cout << "SolverK3BIILU2::GetParameter unsupported operation" << std::endl;
        //throw INMOST::SolverUnsupportedOperation;
        return "";
    }

    void SolverK3BIILU2::SetParameter(std::string name, std::string value) {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
166
167
168
        const char *val = value.c_str();
        if (name == "msglev") solver_data->pParIter->msglev = atoi(val);
        else std::cout << "Parameter " << name << " is unknown" << std::endl;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
    }

    const INMOST_DATA_ENUM_TYPE SolverK3BIILU2::Iterations() const {
        return static_cast<INMOST_DATA_ENUM_TYPE>(SolverIterationNumberK3biilu2(solver_data));
    }

    const INMOST_DATA_REAL_TYPE SolverK3BIILU2::Residual() const {
        return SolverResidualNormK3biilu2(solver_data);
    }

    const std::string SolverK3BIILU2::ReturnReason() const {
        return "Unspecified for K3BIILU2";
    }

    const std::string SolverK3BIILU2::SolverName() const {
        return "k3biilu2";
    }

    void SolverK3BIILU2::Finalize() {
        SolverFinalizeK3biilu2();
    }

    SolverK3BIILU2::~SolverK3BIILU2() {
        this->Clear();
    }
194
195

}