https://www.mediafire.com/file/5qqxes2h ... s.zip/file
gsl library description:
The GNU Scientific Library is a software library for numerical computations in applied mathematics and science. The GSL is written in C; wrappers are available for other programming languages. The GSL is part of the GNU Project and is distributed under the GNU General Public License.
These snippets are starters for those interested.
(I have just started this stuff myself, the .bi files are in the fb distribution already)
code 1
solving linear equations
Code: Select all
#include "gsl/gsl_linalg.bi"
#include "gsl/gsl_blas.bi"
Dim As Double a_data(3,3) = {{ 0.18, 0.60, 0.57, 0.96}, _
{0.41, 0.24, 0.99, 0.58}, _
{0.14, 0.30, 0.97, 0.66}, _
{0.51, 0.13, 0.19, 0.85 }}
Dim As Double copy(3,3) = {{ 0.18, 0.60, 0.57, 0.96}, _
{0.41, 0.24, 0.99, 0.58}, _
{0.14, 0.30, 0.97, 0.66}, _
{0.51, 0.13, 0.19, 0.85 }}
Dim As Double b_data(3) = { 1.0, 2.0, 3.0, 4.0 }
Print "original matrix"
For i As long = 0 To 3
For j As long = 0 To 3
Print copy(i,j);" ";
Next
Print
Next
Print "original right side vector"
For n As Long=0 To 3
Print b_data(n)
Next
print
Dim As gsl_matrix_view m
m = gsl_matrix_view_array (@a_data(0,0), 4, 4)
Dim As gsl_vector_view b
b = gsl_vector_view_array (@b_data(0), 4)
Dim As gsl_vector Ptr x = gsl_vector_alloc (4)
Dim As Long s
Dim As gsl_permutation Ptr p = gsl_permutation_alloc (4)
gsl_linalg_LU_decomp (@m.matrix, p, @s)
gsl_linalg_LU_solve (@m.matrix, p, @b.vector, x)
Print "answer vector"
For n As Long=0 To 3
Print gsl_vector_get (x,n)
Next
gsl_permutation_free (p)
Dim As gsl_matrix Ptr m2 = gsl_matrix_alloc(4,4)
For i As long = 0 To 3
For j As long = 0 To 3
gsl_matrix_set (m2, i, j, copy(i,j))
Next
Next
Dim As gsl_vector Ptr result = gsl_vector_alloc(4)
gsl_blas_dgemv(CblasNoTrans, 1.0,m2,x,0.0,result)
print
Print "check"
printf (!"result = \n")
gsl_vector_fprintf (stdout, result, "%g")
Sleep
matrix multiply
(Just printing the first few.)
Code: Select all
'' Matrix example
'#include "gsl/gsl_matrix.bi"
#include "gsl/gsl_blas.bi"
'' gsl uses the c-style row major order, unlike VB or Fortran
? "A 1024 X 1024 matrix"
Dim As gsl_matrix Ptr m = gsl_matrix_alloc(1024, 1024)
Dim As gsl_matrix Ptr ANS = gsl_matrix_alloc(1024, 1024)
For i As Integer = 0 To 1023
For j As Integer = 0 To 1023
gsl_matrix_set (m, i, j, rnd)
Next
next
For i As Integer = 0 To 2
For j As Integer = 0 To 2
Print "m(";i;",";j;") = "; gsl_matrix_get (m, i, j)
Next
next
?
var t=timer
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m, m, 0.0,ANS)
'gsl_matrix_transpose(m)
print "Please wait . . ."
? "And its mult"
For i As Integer = 0 To 4
For j As Integer = 0 To 4
print gsl_matrix_get (ANS, i, j);" ";
'Print "m(";i;",";j;") = "; gsl_matrix_get (m, i, j);" "
Next
print
next
print timer-t;" seconds"
Sleep
'gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, AT, A, 0.0, ATA);
compare gsl to fb
Code: Select all
#cmdline "-gen gcc -Wc -O2"
#include "gsl/gsl_linalg.bi"
Sub solve1(matrix() As Double, rhs() As Double,solution() As Double)
Dim As Double a_data(Ubound(matrix),Ubound(matrix))
Dim As Double b_data(Ubound(matrix))
For r As Long=0 To Ubound(matrix)
For c As Long=0 To Ubound(matrix)
a_data(r,c)=matrix(r,c)
Next
b_data(r)=rhs(r)
Next
Redim solution(Ubound(matrix))
Dim As gsl_matrix_view m
m = gsl_matrix_view_array (@a_data(0,0),Ubound(matrix)+1, Ubound(matrix)+1)
Dim As gsl_vector_view b
b = gsl_vector_view_array (@b_data(0), Ubound(matrix)+1)
Dim As gsl_vector Ptr x = gsl_vector_alloc (Ubound(matrix)+1)
Dim As Long s
Dim As gsl_permutation Ptr p = gsl_permutation_alloc (Ubound(matrix)+1)
gsl_linalg_LU_decomp (@m.matrix, p, @s)
gsl_linalg_LU_solve (@m.matrix, p, @b.vector, x)
For n As Long=0 To Ubound(matrix)
solution(n)= gsl_vector_get (x,n)
Next
gsl_permutation_free (p)
End Sub
Sub solve2(matrix() As Double, rhs() As Double,solution() As Double)
#Macro setup(All_the_variables)
Dim det As Double=1
Dim mm As Double
Dim m As Long=Ubound(matrix,1)
Dim As Integer x1,k1,g,sign,dimcount
sign=1
Redim l(1 To m,1 To m) As Double
Redim b(1 To m,1 To m) As Double
Redim iv(1 To m) As Double
Redim crhs(1 To m) As Double
Redim solution(1 To m) As Double
#EndMacro
#macro make(copy_matrix_elements)
For i As Long=1 To m
For j As Long=1 To m
b(i,j)=matrix(i,j) 'take a copy of the matrix
If i=j Then
l(i,j)=1.0 'make a unit matrix
Else
l(i,j)=0.0
End If
Next j
crhs(i)=rhs(i)
Next i
#endmacro
#Macro Pivot(If_required)
For x1=1 To Ubound(matrix)
For k1=1 To Ubound(matrix)-1
If k1+1 >Ubound(matrix) Then Exit For
If x1>Ubound(matrix) Then Exit For
If Abs(b(k1+1,1))>Abs(b(x1,1)) Then
sign=-1*sign 'sign changes with each swap
Swap crhs(k1+1),crhs(x1)
For g As Long=1 To Ubound(matrix)
Swap b(k1+1,g),b(x1,g)
Next g
End If
Next k1
Next x1
#EndMacro
#Macro lu(Make_triangular_matrices)
Dim As Long x,y,z
For y=1 To m-1
If b(y,y)=0 Then 'keep ahead of problems by pivoting
pivot(avoid division by zero)
If b(y,y)=0 Then
Print "SINGULAR MATRIX"
Exit Sub
End If
End If
For x=y To m-1
l(x+1,y)=b(x+1,y)/b(y,y) 'l() is lower triangular matrix
mm=l(x+1,y)
b(x+1,y)=0
For z=y+1 To m
b(x+1,z)=b(x+1,z)-mm*b(y,z) 'b() is upper triangular matrix
Next z
Next x
Next y
#EndMacro
#Macro ivector(Make_intermediate_vector)
For n As Long=1 To m
iv(n)=crhs(n)
For j As Long=1 To n-1
iv(n)=iv(n)-l(n,j)*iv(j)
Next j
Next n
#EndMacro
#Macro fvector(Solution_vector)
For n As Long=m To 1 Step -1
solution(n)=iv(n)/b(n,n)
For j As Long = m To n+1 Step -1
solution(n)=solution(n)-(b(n,j)*solution(j)/b(n,n))
Next j
Next n
#EndMacro
' MAIN
setup(All variables) 'set up the working arrays
make(the copy matrix elements)
lu(Find the triangular matrices) 'make upper/lower triangular matrices
ivector(Find the intermediate vector)
fvector(Find the solution)
' END MAIN
End Sub
Randomize 1
Dim Shared As Double mat1(1023,1023)
Dim Shared As Double rhs1(1023)
For r As Long=0 To 1023
For c As Long=0 To 1023
mat1(r,c)=Rnd
Next
rhs1(r)=Rnd
Next
Randomize 1
Dim Shared As Double mat2(1 To 1024,1 To 1024)
Dim Shared As Double rhs2(1 To 1024)
For r As Long=1 To 1024
For c As Long=1 To 1024
mat2(r,c)=Rnd
Next
rhs2(r)=Rnd
Next
Dim As Double t1,t2
Redim As Double answer1(),answer2()
t1=Timer
solve1(mat1(),rhs1(),answer1())
t1=Timer-t1
t2=Timer
solve2(mat2(),rhs2(),answer2())
t2=Timer-t2
For n As Long=0 To 1023
Print answer1(n),answer2(n+1),abs(answer1(n)-answer2(n+1))
Next
Print
Print "gsl ";t1,"mine ";t2
sleep
Some help with blas
https://sites.google.com/a/phys.buruniv ... h-gsl-blas