c++ - 从 C++ 调用 Julia SVD

Question

我尝试使用 Qt C++ 中的 Julia 语言来获取 SVD 函数的结果:

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);

但是当我尝试解析 y 时，我得到了很多垃圾(没关系)，只有 U 和 V 而不是 S:

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

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

输出:

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

那么，我怎样才能从像 SVD 这样的 Julia 函数中正确获取元组呢？

Answer 1

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;
}

编译:

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

运行:

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

输出:

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 输出:

               _
   _       _ _(_)_     |  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>