Search code examples
c++julialinear-algebrasvdijulia-notebook

Call Julia SVD from C++


I try to use Julia language from Qt C++ for getting result of SVD function:

jl_value_t *array_type = jl_apply_array_type(jl_float64_type, 2);
jl_array_t *x  = jl_alloc_array_2d(array_type, matrixForSvd.rows(), matrixForSvd.cols());

jl_value_t * bol = jl_box_bool(true);

double *p = (double*)jl_array_data(x);

int ndims = jl_array_ndims(x);

size_t size0 = jl_array_dim(x,0);
size_t size1 = jl_array_dim(x,1);

// Fill array with data
for(size_t i=0; i<size1; i++)
    for(size_t j=0; j<size0; j++)
        p[j + size0*i] = matrixForSvd(j,i);

jl_function_t *func  = jl_get_function(jl_base_module, "svd");
jl_tupletype_t *y = (jl_tupletype_t*)jl_call1(func, (jl_value_t*)x);

But when I trying to parse y I get many of trash(it's ok) and only U and V but not S:

double *res = (double*)jl_array_data(y);

for(int t = 0; t < 50; t++){
    cout<<res[t]<<endl;
}

Output:

6.94448e-310
1.97626e-323
2.7725e-318
9.88131e-324
9.88131e-324
0
-0.404554 //U
-0.914514
-0.914514
0.404554
0
6.94448e-310
6.94448e-310
1.97626e-323
2.7725e-318
9.88131e-324
9.88131e-324
0
-0.576048 //V
0.817416
-0.817416
-0.576048
0
6.94448e-310
6.94448e-310
1.97626e-323
2.7725e-318

So, how can I correctly get tuple from Julia function like SVD?


Solution

  • svd.cpp:

    #include <julia.h>
    
    #include <iostream>
    #include <vector>
    
    
    void print_vector(jl_array_t *vec) {
        double *data = (double *) jl_array_data(vec);
        size_t n = jl_array_dim(vec, 0);
    
        for(int i = 0; i < n; ++i)
          std::cout << data[i] << " ";
        std::cout << std::endl;
    }
    
    // from column major vector
    void print_2d_matrix(jl_array_t *mat) {
        double *data = (double *) jl_array_data(mat);
        size_t m = jl_array_dim(mat, 0);
        size_t n = jl_array_dim(mat, 1);
    
    
        for(int i = 0; i < m; ++i) {
            for(int j = 0; j < n; ++j)
                std::cout << data[i+j*m] << " ";
            std::cout << std::endl;
        }
    }
    
    int main() {
    
        // matrix represented as a vector in column-major order
        std::vector<double> mat = {
            1, 0, 0, 0,  0, 0, 0, 2,  0, 3, 0, 0,  0, 0, 0, 0,  2, 0, 0, 0
        };
    
        // initialize Julia
        jl_init(NULL);
    
        jl_value_t* dims = (jl_value_t *) jl_eval_string("(4, 5)");
    
        // make sure dims isn't cleaned up by the Julia gc till we're done with it.
        JL_GC_PUSH(&dims);
            // get the svd function
            jl_function_t *svd = jl_get_function(jl_base_module, "svd");
    
            // build a wrapper around the std::vector data to pass our matrix
            // to the svd function
            jl_value_t* array_type = jl_apply_array_type(jl_float64_type, 2);
            jl_array_t *jl_mat = jl_ptr_to_array(array_type, mat.data(), dims, 0);
        JL_GC_POP();
    
        // call svd
        jl_value_t *ret = jl_call1(svd, (jl_value_t*)jl_mat);
    
        // make sure we don't lose our return data
        JL_GC_PUSH1(&ret);
            jl_array_t *jl_U = (jl_array_t*)(jl_fieldref(ret, 0));
            jl_array_t *jl_S = (jl_array_t*)(jl_fieldref(ret, 1));
            jl_array_t *jl_V = (jl_array_t*)(jl_fieldref(ret, 2));
    
            std::cout << "M: " << std::endl;
            print_2d_matrix(jl_mat);
    
            std::cout << "U: " << std::endl;
            print_2d_matrix(jl_U);
    
            std::cout << "S: " << std::endl;
            print_vector(jl_S);
    
            std::cout << "V: " << std::endl;
            print_2d_matrix(jl_V);
    
        JL_GC_POP();
    
        jl_atexit_hook(0);
    
        return 0;
    }
    

    compile with:

    g++ -std=c++14 -fPIC -I$HOME/local/include/julia svd.cpp -L$HOME/local/lib/julia -ljulia
    

    run:

    LD_LIBRARY_PATH=$HOME/local/lib/julia JULIA_HOME=$HOME/local/bin ./a.out
    

    output:

    M: 
    1 0 0 0 2 
    0 0 3 0 0 
    0 0 0 0 0 
    0 2 0 0 0 
    U: 
    0 1 0 0 
    1 0 0 0 
    0 0 0 -1 
    0 0 1 0 
    S: 
    3 2.23607 2 0 
    V: 
    -0 0.447214 -0 0 
    0 0 1 0 
    1 0 0 0 
    -0 0 -0 1 
    0 0.894427 0 0 
    

    julia output:

                   _
       _       _ _(_)_     |  A fresh approach to technical computing
      (_)     | (_) (_)    |  Documentation: http://docs.julialang.org
       _ _   _| |_  __ _   |  Type "?help" for help.
      | | | | | | |/ _` |  |
      | | |_| | | | (_| |  |  Version 0.4.6 (2016-06-19 17:16 UTC)
     _/ |\__'_|_|_|\__'_|  |
    |__/                   |  x86_64-redhat-linux
    
    julia> mat = [1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0]
    4x5 Array{Int64,2}:
     1  0  0  0  2
     0  0  3  0  0
     0  0  0  0  0
     0  2  0  0  0
    
    julia> svd(mat)
    (
    4x4 Array{Float64,2}:
     0.0  1.0  0.0   0.0
     1.0  0.0  0.0   0.0
     0.0  0.0  0.0  -1.0
     0.0  0.0  1.0   0.0,
    
    [3.0,2.23606797749979,2.0,0.0],
    5x4 Array{Float64,2}:
     -0.0  0.447214  -0.0  0.0
      0.0  0.0        1.0  0.0
      1.0  0.0        0.0  0.0
     -0.0  0.0       -0.0  1.0
      0.0  0.894427   0.0  0.0)
    
    julia>