Hi,
I'm trying to use GSL through the FFI, but the program is not linked
to libgslcblas, so it cannot find some functions. Attached a minimal
example.
Cheers, Johan
(use-modules (system foreign)
(rnrs bytevectors))
(define libgsl (dynamic-link "libgsl"))
(define libgslcblas (dynamic-link "libgslcblas"))
;=== Misc. utility functions ========================================
; (make-guard) creates a guard closure that collects objects
; that have memory stored outside the garbage collector with their
; destructors. if the guard is called without arguments, all
; objects in this guard are destroyed.
(define (apply-list l)
(if (null? l)
'()
(begin (apply (caar l) (cdar l))
(apply-list (cdr l)))))
(define (make-guard)
(let ((collection '()))
(case-lambda
(() (apply-list collection))
((destructor obj) (set! collection
(cons (list destructor obj) collection))))))
;====================================================================
;=== matrices =======================================================
(define gsl-matrix-alloc
(pointer->procedure '*
(dynamic-func "gsl_matrix_alloc" libgsl)
(list size_t size_t)))
(define gsl-matrix-free
(pointer->procedure void
(dynamic-func "gsl_matrix_free" libgsl)
(list '*)))
(define (make-gsl-matrix n m guard)
(let ((m (gsl-matrix-alloc n m)))
(guard gsl-matrix-free m)
m))
(define (gsl-matrix->array M)
(let* ((raw (parse-c-struct M (list size_t size_t size_t '* '* int)))
(n (car raw))
(m (cadr raw))
(tda (caddr raw))
(ptr (cadddr raw)))
(make-shared-array
(pointer->bytevector ptr (* tda n) 0 'f64)
(lambda x (list (+ (* tda (car x)) (cadr x))))
n m)))
;====================================================================
;=== vectors ========================================================
(define gsl-vector-alloc
(pointer->procedure '*
(dynamic-func "gsl_vector_alloc" libgsl)
(list size_t)))
(define gsl-vector-free
(pointer->procedure void
(dynamic-func "gsl_vector_free" libgsl)
(list '*)))
(define (gsl-vector->array v)
(let* ((raw (parse-c-struct v (list size_t size_t '* '* int)))
(size (car raw))
(stride (cadr raw))
(ptr (caddr raw)))
(pointer->bytevector ptr size 0 'f64)))
(define (make-gsl-vector n guard)
(let ((v (gsl-vector-alloc n)))
(guard gsl-vector-free v)
v))
;====================================================================
;=== permutations ===================================================
(define permutation-alloc
(pointer->procedure '*
(dynamic-func "gsl_permutation_alloc" libgsl)
(list size_t)))
(define permutation-free
(pointer->procedure void
(dynamic-func "gsl_permutation_free" libgsl)
(list '*)))
(define (make-permutation n guard)
(let ((p (permutation-alloc n)))
(guard permutation-free p)
p))
;====================================================================
;=== linear algebra =================================================
;=== this could be easier, need a c-pointer to a single int ===
(define (dynamic-int)
(make-bytevector (sizeof int)))
(define (dynamic-int->pointer di)
(bytevector->pointer di))
(define (dynamic-int->int di)
(bytevector-sint-ref di 0 (native-endianness) (sizeof int)))
;======
(define linalg-LU-decomp
(pointer->procedure int
(dynamic-func "gsl_linalg_LU_decomp" libgsl)
(list '* '* '*)))
(define linalg-LU-solve
(pointer->procedure int
(dynamic-func "gsl_linalg_LU_solve" libgsl)
(list '* '* '* '*)))
(define linalg-LU-det
(pointer->procedure double
(dynamic-func "gsl_linalg_LU_det" libgsl)
(list '* int)))
(define (make-linalg-det n guard)
(let ((perm (make-permutation n guard))
(sgn (dynamic-int)))
(lambda (m)
(linalg-LU-decomp m perm (dynamic-int->pointer sgn))
(linalg-LU-det m (dynamic-int->int sgn)))))
(define (make-linalg-solve n guard)
(let ((perm (make-permutation n guard))
(sgn (dynamic-int)))
(lambda (A b x)
(linalg-LU-decomp A perm (dynamic-int->pointer sgn))
(linalg-LU-solve A perm b x))))
;====================================================================
(define (randomize-timer)
(let ((time (gettimeofday)))
(set! *random-state*
(seed->random-state (+ (car time)
(cdr time))))))
(let* ((guard (make-guard))
(b (make-gsl-vector 5 guard))
(b-vec (gsl-vector->array b))
(x (make-gsl-vector 5 guard))
(x-vec (gsl-vector->array b))
(A (make-gsl-matrix 5 5 guard))
(A-vec (gsl-matrix->array A))
(det (make-linalg-det 5 guard))
(solve (make-linalg-solve 5 guard)))
(randomize-timer)
(array-index-map! (array-contents A-vec) (lambda X (random:normal)))
(array-index-map! (array-contents b-vec) (lambda X (random:normal)))
(display A-vec) (newline)
(display b-vec) (newline)
(display "=================") (newline)
(display (det A)) (newline)
(solve A b x)
(display x-vec) (newline)
(guard))
