Skip to content
Snippets Groups Projects
Commit a82143c7 authored by liam's avatar liam
Browse files

Add blas1, now working only on gsl-vector (on double floats). Add/fix

unit tests in special-functions.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3049 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 2121c043
No related branches found
No related tags found
No related merge requests found
;********************************************************
; file: blas1.lisp
; description: BLAS level 1
; date: Wed Apr 26 2006 - 15:23
; author: Liam Healy
; modified: Wed Apr 26 2006 - 17:36
;********************************************************
;;; $Id: $
(in-package :gsl)
;;; Currently only includes double-float vector routines.
;;; Not ported: routines that use raw vectors gsl_blas_drotg, gsl_blas_drotmg, gsl_blas_drotm
;;;;****************************************************************************
;;;; Generic
;;;;****************************************************************************
(export '(dot norm asum imax swap copy axpy scal rot))
(defgeneric dot (vector1 vector2)
(:documentation "Dot, or inner, product between vectors."))
(defgeneric norm (x)
(:documentation "The Euclidean norm ||x||_2 = \sqrt {\sum x_i^2} of the vector x."))
(defgeneric asum (x)
(:documentation "The absolute sum \sum |x_i| of the elements of the vector x."))
(defgeneric imax (x)
(:documentation "The index of the largest element of the vector
x. The largest element is determined by its absolute magnitude for
real vectors and by the sum of the magnitudes of the real and
imaginary parts |\Re(x_i)| + |\Im(x_i)| for complex vectors. If the
largest value occurs several times then the index of the first
occurrence is returned."))
(defgeneric swap (x y)
(:documentation "Exchange the elements of the vectors x and y"))
(defgeneric copy (x y)
(:documentation "Copy the elements of the vector x into the vector y."))
(defgeneric axpy (alpha x y)
(:documentation "The sum y = \alpha x + y for the vectors x and y."))
(defgeneric scal (alpha x)
(:documentation "Rescale the vector x by the multiplicative factor alpha."))
(defgeneric rot (x y c s)
(:documentation
"Apply a Givens rotation (x', y') = (c x + s y, -s x + c y) to the vectors x, y."))
;;;;****************************************************************************
;;;; Double vectors
;;;;****************************************************************************
(defun-gsl dot (((pointer vec1) :pointer) ((pointer vec2) :pointer))
"gsl_blas_ddot"
:method
((vec1 gsl-vector) (vec2 gsl-vector))
:return (:double))
(defun-gsl norm (((pointer vec) :pointer))
"gsl_blas_dnrm2"
:method
((vec gsl-vector))
:c-return-value :return
:return (:double))
(defun-gsl asum (((pointer vec) :pointer))
"gsl_blas_dasum"
:method
((vec gsl-vector))
:c-return-value :return
:return (:double))
(defun-gsl imax (((pointer vec) :pointer))
"gsl_blas_idamax"
:method ((vec gsl-vector))
:c-return-value :return
:return (:int))
(defun-gsl swap (((pointer vec1) :pointer) ((pointer vec2) :pointer))
"gsl_blas_dswap"
:method ((vec1 gsl-vector) (vec2 gsl-vector))
:after ((cl-invalidate vec2)))
(defun-gsl copy (((pointer vec1) :pointer) ((pointer vec2) :pointer))
"gsl_blas_dcopy"
:method ((vec1 gsl-vector) (vec2 gsl-vector))
:after ((cl-invalidate vec2)))
(defun-gsl axpy
((alpha :double) ((pointer vec1) :pointer) ((pointer vec2) :pointer))
"gsl_blas_daxpy"
:method (alpha (vec1 gsl-vector) (vec2 gsl-vector))
:after ((cl-invalidate vec2)))
(defun-gsl scal ((alpha :double) ((pointer vec) :pointer))
"gsl_blas_dscal"
:method (alpha (vec gsl-vector))
:after ((cl-invalidate vec))
:c-return-value :void)
(defun-gsl rot
(((pointer vec1) :pointer) ((pointer vec2) :pointer) (c :double) (s :double))
"gsl_blas_drot"
:method ((vec1 gsl-vector) (vec2 gsl-vector) (c double-float) (s double-float))
:after ((cl-invalidate vec1 vec2)))
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
......@@ -3,7 +3,7 @@
; description: Using GSL storage.
; date: Sun Mar 26 2006 - 16:32
; author: Liam M. Healy
; modified: Wed Apr 19 2006 - 00:46
; modified: Wed Apr 26 2006 - 17:37
;********************************************************
;;; $Id: $
......@@ -184,6 +184,7 @@
;;;;****************************************************************************
(defun cl-invalidate (&rest objects)
"Mark the CL image of the GSL array as invalid."
(mapc (lambda (obj) (setf (cl-invalid obj) t))
objects))
......
......@@ -3,7 +3,7 @@
; description: Definition of GSLL system
; date:
; author: Liam Healy
; modified: Tue Apr 25 2006 - 14:56
; modified: Wed Apr 26 2006 - 16:20
;********************************************************
;;; $Id: $
......@@ -19,6 +19,7 @@
:components
((:file "init")
(:file "interface" :depends-on (init))
;; http://www.cs.northwestern.edu/academics/courses/325/readings/lisp-unit.html
(:file "lisp-unit")
(:file "tests")))
(:module general
......@@ -37,7 +38,8 @@
(:file "vector" :depends-on (data))
(:file "matrix" :depends-on (data))
(:file "permutation" :depends-on (data))
(:file "combination" :depends-on (data))))
(:file "combination" :depends-on (data))
(:file "blas1")))
(:module special-functions
:depends-on (init)
:components
......
......@@ -3,7 +3,7 @@
; description: Exponential integrals
; date: Tue Mar 21 2006 - 17:37
; author: Liam M. Healy
; modified: Sat Mar 25 2006 - 22:11
; modified: Wed Apr 26 2006 - 10:13
;********************************************************
;;; $Id: $
......@@ -86,3 +86,34 @@
:documentation
"The Arctangent integral, which is defined as @math{AtanInt(x) = \int_0^x dt \arctan(t)/t}."
:return (sf-result))
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
(lisp-unit:define-test exponential-integrals
(lisp-unit:assert-error 'gsl-error (expint-E1 0.0d0))
(lisp-unit:assert-first-fp-equal
"0.219383934396d+00"
(expint-E1 1.0d0))
(lisp-unit:assert-first-fp-equal
"0.495423435600d+01"
(expint-Ei 2.0d0))
(lisp-unit:assert-first-fp-equal
"0.136373067344d+01"
(Shi 1.25d0))
(lisp-unit:assert-first-fp-equal
"0.121731730091d+01"
(Chi 1.25d0))
(lisp-unit:assert-first-fp-equal
"0.868892654126d+00"
(expint-3 1.25d0))
(lisp-unit:assert-first-fp-equal
"0.114644641567d+01"
(si 1.25d0))
(lisp-unit:assert-first-fp-equal
"0.434300724034d+00"
(ci 1.25d0))
(lisp-unit:assert-first-fp-equal
"0.110361916168d+01"
(atanint 1.25d0)))
......@@ -3,7 +3,7 @@
; description: Fermi-Dirac function.
; date: Sat Apr 22 2006 - 16:12
; author: Liam M. Healy
; modified: Sat Apr 22 2006 - 17:41
; modified: Wed Apr 26 2006 - 09:51
;********************************************************
;;; $Id: $
......@@ -103,4 +103,4 @@
(fermi-dirac-1/2 2.0d0))
(lisp-unit:assert-first-fp-equal
"0.170141327798d+01"
(fp-string (fermi-dirac-inc-0 2.0d0 0.5d0))))
(fermi-dirac-inc-0 2.0d0 0.5d0)))
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment