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

namespace INMOST {

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

Dmitry Bagaev's avatar
Dmitry Bagaev committed
12
    SolverInterface *SolverK3BIILU2::Copy(const SolverInterface *other) {
13 14 15 16 17 18 19 20 21
        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
22 23 24 25 26
        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
27
        return this;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
28 29 30
    }

    void SolverK3BIILU2::Assign(const SolverInterface *other) {
31 32 33 34 35 36 37 38 39
        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
40 41 42 43 44 45 46 47 48 49 50
        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
51 52
    void SolverK3BIILU2::Setup(int *argc, char ***argv, SolverParameters &p) {
        SolverInitDataK3biilu2(&solver_data, communicator, p.solverPrefix.c_str());
53 54 55 56 57 58 59
        if (p.internalFile != "") {
            SolverInitializeK3biilu2(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);
            }
        }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
    }

    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;
76
#if defined(USE_MPI)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
77 78
            MPI_Comm_size(A.GetCommunicator(), &nproc);
            MPI_Comm_rank(A.GetCommunicator(), &myid);
79
#else
Dmitry Bagaev's avatar
Dmitry Bagaev committed
80 81
            nproc = 1;
            myid = 0;
82
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
83 84 85
            int *ibl = (int *) malloc(sizeof(int) * (nproc + 1));
            ibl[0] = 0;
            int n = A.Size();
86
#if defined(USE_MPI)
87
            INMOST_DATA_ENUM_TYPE mbeg, mend;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
88 89 90 91
            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());
92
#else
Dmitry Bagaev's avatar
Dmitry Bagaev committed
93
            ibl[1] = n;
94
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
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 128 129 130 131 132
            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;
133
        RHS.GetInterval(vbeg, vend);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
134 135
        vector_k3biilu2 *rhs_data = NULL;
        vector_k3biilu2 *solution_data = NULL;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
136 137 138 139 140 141 142 143 144 145 146 147 148 149 150

        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]);
151
        iter_time = solver_data->dstat[9];
152
        time_prec = solver_data->dstat[7];
Dmitry Bagaev's avatar
Dmitry Bagaev committed
153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
        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 {
169 170 171 172 173
        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]);
174
#if !defined(SILENCE_SET_PARAMETER)
175
        std::cout << "Parameter " << name << " is unknown" << std::endl;
176
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
177 178 179 180
        return "";
    }

    void SolverK3BIILU2::SetParameter(std::string name, std::string value) {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
181 182
        const char *val = value.c_str();
        if (name == "msglev") solver_data->pParIter->msglev = atoi(val);
183 184 185 186 187 188
        else if (name == "kovl") solver_data->pParams->ncycle = atoi(val);
        else if (name == "tau") {
            solver_data->pParams->tau1 = atof(val);
            solver_data->pParams->tau2 = -1.0;
        } else if (name == "nit") solver_data->pParIter->maxit = atoi(val);
        else if (name == "eps") solver_data->pParIter->eps = atof(val);
189
#if !defined(SILENCE_SET_PARAMETER)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
190
        else std::cout << "Parameter " << name << " is unknown" << std::endl;
191
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
192 193
    }

194
    INMOST_DATA_ENUM_TYPE SolverK3BIILU2::Iterations() const {
195
        return static_cast<INMOST_DATA_ENUM_TYPE>(solver_data->istat[2]);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
196 197
    }

198
    INMOST_DATA_REAL_TYPE SolverK3BIILU2::Residual() const {
199
        return solver_data->dstat[2];
Dmitry Bagaev's avatar
Dmitry Bagaev committed
200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216
    }

    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();
    }
217

218
}