RosettaCodeData/Task/Singular-value-decomposition/FreeBASIC/singular-value-decompositio...

65 lines
1.3 KiB
Plaintext

#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