#include once "inc\gsl\gsl_linalg.bi"

Sub MatrixPrint(r As Integer, c As Integer, m() As Double)
    For i As Integer = 0 To r - 1
        Print "|";
        For j As Integer = 0 To c - 1
            Print Using "#############.########## "; m(i * c + j);
        Next j
        Print Chr(8); "|"
    Next i
    Print
End Sub

Dim As Integer r = 2
Dim As Integer c = 2
Dim As Integer l = r * c

Dim As Double a(l - 1)
a(0) = 3: a(1) = 0: a(2) = 4: a(3) = 5
Dim As Double v(l - 1)
v(0) = 0: v(1) = 0: v(2) = 0: v(3) = 0
Dim As Double s(c - 1)
s(0) = 0: s(1) = 0

Dim As gsl_matrix Ptr a_mat = gsl_matrix_alloc(r, c)
Dim As gsl_matrix Ptr v_mat = gsl_matrix_alloc(c, c)
Dim As gsl_vector Ptr s_vec = gsl_vector_alloc(c)
Dim As gsl_vector Ptr work = gsl_vector_alloc(c)

Dim As Integer i, j
For i = 0 To r - 1
    For j = 0 To c - 1
        gsl_matrix_set(a_mat, i, j, a(i * c + j))
    Next j
Next i

gsl_linalg_SV_decomp(a_mat, v_mat, s_vec, work)

For i = 0 To r - 1
    For j = 0 To c - 1
        a(i * c + j) = gsl_matrix_get(a_mat, i, j)
        v(i * c + j) = gsl_matrix_get(v_mat, i, j)
    Next j
Next i

For i = 0 To c - 1
    s(i) = gsl_vector_get(s_vec, i)
Next i

Print "U:"
MatrixPrint(r, c, a())

Print "S:"
MatrixPrint(1, c, s())

Print "VT:"
MatrixPrint(r, c, v())

gsl_matrix_free(a_mat)
gsl_matrix_free(v_mat)
gsl_vector_free(s_vec)
gsl_vector_free(work)

Sleep