sparse_complex 0.1.4

A simple solver for complex sparse matrices.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
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

#include <complex>
#include <Eigen/Sparse>
#include <Eigen/SparseLU>
#include <vector>
#include <iostream>

typedef std::complex<double> cdouble;
typedef std::complex<float> cfloat;
typedef Eigen::Triplet<cdouble> T;
typedef Eigen::Triplet<cfloat> T32;

extern "C" struct cmplx64
{
    double re;
    double im;
};

extern "C" struct cmplx32
{
    float re;
    float im;
};

extern "C" void solve_cpp(const cmplx64 *values, const size_t *rows, const size_t *cols, const size_t n_value, cmplx64 *b, const size_t size)
{
    Eigen::SparseMatrix<cdouble> A(size, size);
    A.reserve(n_value);

    std::vector<T> triplets;
    triplets.reserve(n_value);

    Eigen::VectorXcd b_(size);
    Eigen::VectorXcd x(size);

    Eigen::SparseLU< Eigen::SparseMatrix<cdouble> > solver;

    double re, im;
    cdouble value;
    size_t row, col, i;

    for (i = 0; i < n_value; i++)
    {
        re = (values + i)->re;
        im = (values + i)->im;
        value = cdouble(re, im);
        row = *(rows + i);
        col = *(cols + i);
        triplets.push_back(T(row, col, value));
    }

    for (i = 0; i < size; i++)
    {

        re = (b + i)->re;
        im = (b + i)->im;
        value = cdouble(re, im);
        b_[i] = value;
    }

    A.setFromTriplets(triplets.begin(), triplets.end());

    solver.compute(A);
    x = solver.solve(b_);

    for (size_t i = 0; i < size; i++)
    {
        re = x[i].real();
        im = x[i].imag();
        (b + i)->re = re;
        (b + i)->im = im;
    }
}

extern "C" void solve_cpp32(const cmplx32 *values, const size_t *rows, const size_t *cols, const size_t n_value, cmplx32 *b, const size_t size)
{
    Eigen::SparseMatrix<cfloat> A(size, size);
    A.reserve(n_value);

    std::vector<T32> triplets;
    triplets.reserve(n_value);

    Eigen::VectorXcf b_(size);
    Eigen::VectorXcf x(size);

    Eigen::SparseLU< Eigen::SparseMatrix<cfloat> > solver;

    float re, im;
    cfloat value;
    size_t row, col, i;

    for (i = 0; i < n_value; i++)
    {
        re = (values + i)->re;
        im = (values + i)->im;
        value = cfloat(re, im);
        row = *(rows + i);
        col = *(cols + i);
        triplets.push_back(T32(row, col, value));
    }

    for (i = 0; i < size; i++)
    {

        re = (b + i)->re;
        im = (b + i)->im;
        value = cfloat(re, im);
        b_[i] = value;
    }

    A.setFromTriplets(triplets.begin(), triplets.end());

    solver.compute(A);
    x = solver.solve(b_);

    for (size_t i = 0; i < size; i++)
    {
        re = x[i].real();
        im = x[i].imag();
        (b + i)->re = re;
        (b + i)->im = im;
    }
}