SolverSUPERLU.cpp 5.27 KB
Newer Older
1
2
3
4
#include "SolverSUPERLU.h"

namespace INMOST {

Dmitry Bagaev's avatar
Dmitry Bagaev committed
5
6
7
8
9
10
11
12
    SolverSUPERLU::SolverSUPERLU() {
        set_default_options(&options);
        StatInit(&stat);
        perm_c = NULL;
        perm_r = NULL;
        a_size = 0;
    }

Dmitry Bagaev's avatar
Dmitry Bagaev committed
13
    SolverInterface *SolverSUPERLU::Copy(const SolverInterface *other) {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
14
15
16
17
18
19
20
        throw INMOST::SolverUnsupportedOperation;
    }

    void SolverSUPERLU::Assign(const SolverInterface *other) {
        throw INMOST::SolverUnsupportedOperation;
    }

Dmitry Bagaev's avatar
Dmitry Bagaev committed
21
    void SolverSUPERLU::Setup(int *argc, char ***argv, SolverParameters &p) {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
        //read options from file and arguments
    }

    void SolverSUPERLU::SetMatrix(Sparse::Matrix &A, bool ModifiedPattern, bool OldPreconditioner) {
        //check that the run is serial!
        int *ia, *ja, nnz = 0;
        double *a;
        int mbeg = A.GetFirstIndex();
        int mend = A.GetLastIndex();
        int size = 0;
        remap = new int[mend - mbeg];
        for (int k = 0; k < mend - mbeg; ++k) {
            Sparse::Row &r = A[k + mbeg];
            if (r.Size()) {
                double nrm = 0;
                for (INMOST_DATA_ENUM_TYPE l = 0; l < r.Size(); ++l) {
                    nrm += r.GetValue(l) * r.GetValue(l);
                }
                if (nrm) {
                    std::sort(r.Begin(), r.End());
                    nnz += r.Size();
                    remap[k] = size;
                    size++;
                } else {
                    remap[k] = -1;
                }
            } else {
                remap[k] = -1;
            }
51
        }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
        ia = (int *) malloc(sizeof(int) * (size + 1));
        ja = (int *) malloc(sizeof(int) * nnz);
        a = (double *) malloc(sizeof(double) * nnz);
        int q = 0, f = 0;
        ia[0] = 0;
        for (int k = 0; k < mend - mbeg; ++k) {
            if (remap[k] != -1) {
                Sparse::Row &r = A[k + mbeg];
                for (INMOST_DATA_ENUM_TYPE l = 0; l < r.Size(); ++l) {
                    if (remap[r.GetIndex(l) - mbeg] != -1) {
                        ja[q] = remap[r.GetIndex(l) - mbeg];
                        a[q] = r.GetValue(l);
                        ++q;
                    } else { //if( fabs(a[q]) > 1.0e-9 )
                        std::cout << "Matrix has connections to singular rows" << std::endl;
                    }
                }
                ia[f + 1] = q;
                f++;
            }
72
        }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
73
74
75
76
77
78
79
80
81
82
83
84
85
        dCreate_CompCol_Matrix(&(this->A), size, size, nnz, a, ja, ia, SLU_NR, SLU_D, SLU_GE);
        a_size = size;
        perm_c = new int[size];
        perm_r = new int[size];
    }

    bool SolverSUPERLU::Solve(INMOST::Sparse::Vector &RHS, INMOST::Sparse::Vector &SOL) {
        double *inout = new double[a_size];
        int mbeg = RHS.GetFirstIndex(), mend = RHS.GetLastIndex();
        for (int k = 0; k < mend - mbeg; ++k) {
            if (remap[k] != -1) {
                inout[remap[k]] = RHS[k + mbeg];
            }
86
        }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
87
88
89
90
91
92
93
94
95
        SuperMatrix B;
        dCreate_Dense_Matrix(&B, a_size, 1, inout, a_size, SLU_DN, SLU_D, SLU_GE);
        dgssv(&options, &A, perm_c, perm_r, &L, &U, &B, &stat, &info);
        Destroy_SuperMatrix_Store(&B);
        bool ret = (info == 0);
        for (int k = 0; k < mend - mbeg; ++k) {
            if (remap[k] != -1) {
                SOL[k + mbeg] = inout[remap[k]];
            }
96
        }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
        delete[] inout;
        return ret;
    }

    bool SolverSUPERLU::Clear() {
        Destroy_CompCol_Matrix(&A);
        Destroy_SuperNode_Matrix(&L);
        Destroy_CompCol_Matrix(&U);
        StatFree(&stat);
        if (perm_c != NULL) delete[] perm_c;
        if (perm_r != NULL) delete[] perm_r;
        if (remap != NULL) delete[] remap; //allocated outside
        return true;
    }

    bool SolverSUPERLU::isMatrixSet() {
        return a_size != 0;
    }

    std::string SolverSUPERLU::GetParameter(std::string name) const {
Kirill Terekhov's avatar
Kirill Terekhov committed
117
#if !defined(SILENCE_SET_PARAMETER)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
118
        std::cout << "SolverSUPERLU::GetParameter unsupported operation" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
119
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
120
121
122
123
124
        //throw INMOST::SolverUnsupportedOperation;
        return "";
    }

    void SolverSUPERLU::SetParameter(std::string name, std::string value) {
Kirill Terekhov's avatar
Kirill Terekhov committed
125
#if !defined(SILENCE_SET_PARAMETER)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
126
        std::cout << "SolverSUPERLU::SetParameter unsupported operation" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
127
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
128
129
130
        //throw INMOST::SolverUnsupportedOperation;
    }

131
    INMOST_DATA_ENUM_TYPE SolverSUPERLU::Iterations() const {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
132
133
134
        return 1;
    }

135
    INMOST_DATA_REAL_TYPE SolverSUPERLU::Residual() const {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
136
137
138
139
140
141
142
143
144
145
146
147
148
        return 0;
    }

    const std::string SolverSUPERLU::ReturnReason() const {
        char reason_str[256];
        if (info <= a_size) {
            sprintf(reason_str, "diagonal entry of U-factor is exactly singular at %d/%d", info, a_size);
        } else if (info > a_size) {
            sprintf(reason_str, "memory allocation failed after %d bytes were allocated", info - a_size);
        } else if (info == 0) {
            strcpy(reason_str, "factorization was successfull");
        } else {
            sprintf(reason_str, "unknown exit code %d", info);
149
        }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
150
151
        return std::string(reason_str);
    }
152

Dmitry Bagaev's avatar
Dmitry Bagaev committed
153
154
155
    const std::string SolverSUPERLU::SolverName() const {
        return "superlu";
    }
156

Dmitry Bagaev's avatar
Dmitry Bagaev committed
157
158
159
    void SolverSUPERLU::Finalize() {
        //nothing to do here
    }
160

Dmitry Bagaev's avatar
Dmitry Bagaev committed
161
    SolverSUPERLU::~SolverSUPERLU() {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
162
        //this->Clear();
Dmitry Bagaev's avatar
Dmitry Bagaev committed
163
    }
164
165
166


}