Using gsl library

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
Post Reply
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Using gsl library

Post by dodicat »

The 64 bit dll files
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
  
code 2
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);
  
code 3
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

 
Please put the dll and dll.a files into the same folder as these tests, or vice versa.
Some help with blas
https://sites.google.com/a/phys.buruniv ... h-gsl-blas
Avata
Posts: 102
Joined: Jan 17, 2021 7:27

Re: Using gsl library

Post by Avata »

Where could download the static library file(.dll and .a)
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Using gsl library

Post by dodicat »

try srvaldez's link
viewtopic.php?t=28508
I have only used the dll so far, let me know if the static works please.
Post Reply