Skip to content
Snippets Groups Projects
Commit 4647599b authored by Liam Healy's avatar Liam Healy
Browse files

Pass complex scalars to GSL functions

By sending the two parts of a complex as separate doubles (complex
double-float), or packing the two parts into a one double-float
argument (complex single-float), it is possible to pass complex
numbers in GSL function calls.  This is very likely to be
non-portable, so a test is conducted and the results
in *pass-complex-scalar-as-two-reals* determine whether defmfun does
this translation or just creates a function that signals an error.
The complex tests have been modified so that no longer expect an
error:
SBCL 64: TOTAL: 1223 assertions passed, 5 failed, 0 execution errors.
CCL 64: TOTAL: 1228 assertions passed, 0 failed, 0 execution errors.
This also works in 32 bit SBCL.
parent d8fd822a
No related branches found
No related tags found
No related merge requests found
;; Functions of complex numbers
;; Liam Healy 2009-01-13 21:19:38EST complex.lisp
;; Time-stamp: <2009-01-13 21:57:29EST complex.lisp>
;; $Id: $
(in-package :gsl)
;;;;****************************************************************************
;;; Properties
;;;;****************************************************************************
(defmfun argument (number)
"gsl_complex_arg"
((number complex-double-c))
:c-return :double
:documentation
"The angle of the complex number.")
(defmfun modulus (number)
"gsl_complex_abs"
((number complex-double-c))
:c-return :double
:documentation "The magnitude, or modulus of the complex number.")
(defmfun abs2 (number)
"gsl_complex_abs2"
((number complex-double-c))
:c-return :double
:documentation "The magnitude squared of the complex number.")
(defmfun logabs2 (number)
"gsl_complex_logabs2"
((number complex-double-c))
:c-return :double
:documentation "The magnitude squared of the complex number.")
;;;;****************************************************************************
;;; Complex arithmetic
;;;;****************************************************************************
;;; The GSL functions provided duplicate what CL has built-in.
;; Comparison of floating point numbers using sequence integers
;; Liam Healy 2008-01-22 19:00:17EST floating-point.lisp
;; Time-stamp: <2008-04-04 11:35:20EDT floating-point.lisp>
;; Time-stamp: <2009-01-14 21:33:50EST floating-point.lisp>
;; $Id$
;;; All floats can be represented by integers. There are two slightly
......@@ -76,7 +76,9 @@
(declare (ignore expbits))
(when (eq ieee754 t) (setf ieee754 (type-of float)))
(if ieee754
(* (dpb 1 (ieee754-sign-bit ieee754) (+ (ash exp sigbits) mant)))
(dpb (if (minusp sign) 1 0)
(ieee754-sign-bit ieee754)
(+ (ash exp sigbits) mant))
(* sign (+ (ash exp sigbits) mant)))))
;; (float-as-integer most-negative-single-float)
......
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2009-01-07 22:09:19EST gsll.asd>
;; Time-stamp: <2009-01-14 22:12:10EST gsll.asd>
;; $Id$
(asdf:defsystem "gsll"
......@@ -17,7 +17,8 @@
(:file "conditions" :depends-on (init))
(:file "mobject" :depends-on (init))
(:file "types" :depends-on (init))
(:file "element-types" :depends-on (init types))
(:file "complex-types" :depends-on (types))
(:file "element-types" :depends-on (init complex-types))
(:file "number-conversion" :depends-on (init))
(:file "interface"
:depends-on (init conditions element-types number-conversion))
......@@ -29,7 +30,8 @@
(:module floating-point
:depends-on (init)
:components
((:file "ieee-modes")))
((:file "ieee-modes")
(:file "floating-point")))
(:file "mathematical" :depends-on (init))
;; complex numbers not necessary? Just make a struct.
(:module data
......
;; Complex number types
;; Liam Healy 2009-01-13 21:24:05EST complex-types.lisp
;; Time-stamp: <2009-01-14 22:45:53EST complex-types.lisp>
;; $Id: $
(in-package :gsl)
;;;;****************************************************************************
;;;; Complex types
;;;;****************************************************************************
;;; GSL represents complex variables as a struct of two reals and
;;; passes them by value (not a pointer) to functions. Few FFIs have
;;; the capability of handling that. However, it does seem to work
;;; (SBCL and CCL at least) that we can treat each complex argument as
;;; a succession of two double-float arguments, or two single-floats
;;; packed into a double float. Here we probe that the first is true;
;;; if so, functions will automatically be converted.
(defparameter *pass-complex-scalar-as-two-reals*
(eql 5.0d0
(ignore-errors
(cffi:foreign-funcall "gsl_complex_abs" :double 3.0d0 :double 4.0d0 :double)))
"A complex number can be passed as two adjacent reals.")
;;; GSL defines complex numbers in a struct, and passes the struct by
;;; value. CFFI does not support call by value for structs, so we
;;; cannot use functions that call or return complex scalars.
;;; See /usr/include/gsl/gsl_complex.h
(cffi:defcstruct complex-float-c
(dat :float :count 2))
(cffi:defcstruct complex-double-c
(dat :double :count 2))
#+long-double
(cffi:defcstruct complex-long-double-c
(dat :long-double :count 2))
(defun clean-type (type)
;; SBCL at least will specify limits on the type, e.g.
;; (type-of #C(1.0 2.0))
;; (COMPLEX (DOUBLE-FLOAT 1.0 2.0))
;; This cleans that up to make
;; (clean-type (type-of #C(1.0 2.0)))
;; (COMPLEX DOUBLE-FLOAT)
(if (and (subtypep type 'complex) (listp (second type)))
(list (first type) (first (second type)))
type))
(defun component-float-type (eltype)
"The type of the component of this type (complex)."
(if (subtypep eltype 'complex)
;; complex: use the component type
(second eltype)
eltype))
(define-condition pass-complex-by-value (error)
()
(:report
(lambda (condition stream)
(declare (ignore condition))
(format stream "Cannot pass complex scalars to GSL functions.")))
(:documentation
"An error indicating that this implementation and platform are
unable to pass complex numbers to the GSL libary by value."))
(defun pack-complex-as-double (number)
"Pack a number of type (complex single-float) into a double-float,
which may be accepted as a struct of type gsl_complex_float.
This is hideously non-portable."
(integer-as-float
(dpb
(float-as-integer (imagpart number) t)
(byte 32 32)
(dpb
(float-as-integer (realpart number) t)
(byte 32 0)
0))
'double-float))
(defun passing-complex-by-value (cfind)
"Substitution in defmfun so that complex numbers can be passed by
value to the GSL structs defined above. This is non-portable."
(if (eq (component-float-type (cffi-cl (third cfind))) 'single-float)
;; single-float
`(:double (pack-complex-as-double ,(first cfind)))
;; double-float
`(:double
(realpart ,(first cfind))
:double
(imagpart ,(first cfind)))))
;; Helpers that define a single GSL function interface
;; Liam Healy 2009-01-07 22:02:20EST defmfun-single.lisp
;; Time-stamp: <2009-01-11 22:37:31EST defmfun-single.lisp>
;; Time-stamp: <2009-01-14 22:03:12EST defmfun-single.lisp>
;; $Id: $
(in-package :gsl)
......@@ -148,52 +148,50 @@
outputs
(unless (eq c-return :void)
(list cret-name)))))
(wrap-letlike
(or allocated-decl complex-args)
(mapcar #'wfo-declare
(append allocated-decl
(mapcar #'rest complex-args)))
'cffi:with-foreign-objects
`(#-native
,@(mapcar (lambda (v) `(copy-cl-to-c ,v)) inputs)
,@before
(let ((,cret-name
(cffi:foreign-funcall
,gsl-name
,@(mapcan
(lambda (arg)
(let ((cfind ; variable is complex
(first (member (st-symbol arg)
complex-args :key 'first))))
(if cfind ; so substitute call to complex-to-gsl
`(,(third cfind)
(complex-to-gsl ,(first cfind) ,(second cfind)))
;; otherwise use without conversion
(list (if (member (st-symbol arg) allocated)
:pointer
(st-type arg))
(st-symbol arg)))))
c-arguments)
,cret-type)))
,@(case c-return
(:void `((declare (ignore ,cret-name))))
(:error-code ; fill in arguments
`((check-gsl-status ,cret-name
',(or (defgeneric-method-p name) name)))))
#-native
,@(when outputs
(mapcar
(lambda (x) `(setf (cl-invalid ,x) t (c-invalid ,x) nil))
outputs))
,@(when (eq cret-type :pointer)
`((check-null-pointer
,cret-name
,@'('memory-allocation-failure "No memory allocated."))))
,@after
(values
,@(defmfun-return
c-return cret-name clret allocated return return-supplied-p
enumeration outputs))))))))
(if (and complex-args (not *pass-complex-scalar-as-two-reals*))
'(error 'pass-complex-by-value) ; arglist should be declared ignore
(wrap-letlike
allocated-decl
(mapcar #'wfo-declare allocated-decl)
'cffi:with-foreign-objects
`(#-native
,@(mapcar (lambda (v) `(copy-cl-to-c ,v)) inputs)
,@before
(let ((,cret-name
(cffi:foreign-funcall
,gsl-name
,@(mapcan
(lambda (arg)
(let ((cfind ; variable is complex
(find (st-symbol arg) complex-args :key 'first)))
(if cfind ; make two successive scalars
(passing-complex-by-value cfind)
;; otherwise use without conversion
(list (if (member (st-symbol arg) allocated)
:pointer
(st-type arg))
(st-symbol arg)))))
c-arguments)
,cret-type)))
,@(case c-return
(:void `((declare (ignore ,cret-name))))
(:error-code ; fill in arguments
`((check-gsl-status ,cret-name
',(or (defgeneric-method-p name) name)))))
#-native
,@(when outputs
(mapcar
(lambda (x) `(setf (cl-invalid ,x) t (c-invalid ,x) nil))
outputs))
,@(when (eq cret-type :pointer)
`((check-null-pointer
,cret-name
,@'('memory-allocation-failure "No memory allocated."))))
,@after
(values
,@(defmfun-return
c-return cret-name clret allocated return return-supplied-p
enumeration outputs)))))))))
(defun defmfun-return
(c-return cret-name clret allocated return return-supplied-p enumeration outputs)
......
;; Mapping of element type names
;; Liam Healy 2008-04-13 11:22:46EDT element-types.lisp
;; Time-stamp: <2008-12-31 21:10:07EST element-types.lisp>
;; Time-stamp: <2009-01-13 21:31:12EST element-types.lisp>
;; $Id$
;;; The different element type forms:
;;; C standard full, or "CFFI" :unsigned-char
;;; GSL splice name "uchar"
;;; C explicit :uint8 (not used here)
;;; CL '(unsigned-byte 8)
;;; Single 'unsigned-byte-8
;;; Functions to perform conversions
;;; CL -> single in function #'cl-single
;;; CL -> GSL in function #'cl-gsl
;;; CL -> CFFI (Cstd) in function #'cl-cffi
;;; CFFI (Cstd) -> CL in function #'cffi-cl
;;; Sources of equivalence
;;; Cstd -> CL in alist *cstd-cl-type-mapping*
;;; Cstd -> GSL in alist *cstd-gsl-mapping*
;;; :long-double may be pushed onto *features* if
;;; the implementation supports long doubles in CFFI.
(in-package :gsl)
;;;;****************************************************************************
;;;; Basic definition
;;;;****************************************************************************
;;; Preliminary definitions. These are not used directly outside this
;;; file; they exist to define the two variables
;;; *cstd-cl-type-mapping* and *cstd-gsl-mapping* which are used by
;;; the conversion functions.
(defparameter *cstd-integer-types*
'(:char :unsigned-char
:short :unsigned-short
:int :unsigned-int
:long :unsigned-long
#+int64
:long-long
#+int64
:unsigned-long-long)
;; http://common-lisp.net/project/cffi/manual/html_node/Built_002dIn-Types.html
"List of integer types supported by CFFI, from the CFFI docs.")
(defparameter *fp-type-mapping*
'((:float . single-float) (:double . double-float)
;; For those implementations that support a separate long-double
;; type assume this mapping:
#+long-double (:long-double . long-float)
(complex-float-c . (complex single-float))
(complex-double-c . (complex double-float))
;; For those implementations that support a separate long-double
;; type assume this mapping:
#+long-double
(complex-long-double-c . (complex long-float)))
;; Ordered by: real shortest to longest, then complex shortest to longest.
"List of floating point types supported by CFFI from the CFFI docs
plus corresponding complex types.")
(defparameter *gsl-splice-int-types*
;; list | grep -i 'gsl_vector.*_alloc\b'
'("char" "int" "long" "short" "uchar" "uint" "ulong" "ushort")
"The list of integer types that can be spliced into function names.")
(defparameter *gsl-splice-fp-types*
;; list | grep -i 'gsl_vector.*_alloc\b'
;; Ordered by: real shortest to longest, then complex shortest to longest.
;; Must match *fp-type-mapping*.
'("float" "" #+long-double "long_double"
"complex_float" "complex" #+long-double "complex_long_double")
"The list of floating point types that can be spliced into function names.")
;;; Mapping alists used by conversion functions.
;;; Used by #'cl-gsl
(defparameter *cstd-cl-type-mapping*
(append
*fp-type-mapping*
(mapcar
(lambda (type)
(cons
type
(list
(if (string-equal type "uns" :end1 3)
'unsigned-byte
'signed-byte)
(* 8 (cffi:foreign-type-size type)))))
*cstd-integer-types*))
;; Be careful when reverse associating, as there may be several C
;; types that map to a single CL type.
"An alist of the C standard types as keywords, and the CL type
The exception is complex types, which don't have a definition
in the C standard; in that case, the C type is the GSL struct
definition.")
(defmacro floating-point-association (splice-list)
`(mapcar
#'cons
(mapcar #'first *fp-type-mapping*)
(subseq ,splice-list 0 (length *fp-type-mapping*))))
(defparameter *cstd-gsl-mapping*
(append
;; The integer types
(remove-if-not
(lambda (x) (find (rest x) *gsl-splice-int-types* :test 'string-equal))
(mapcar
(lambda (type)
(cons type
(let ((ut
(if (string-equal type "uns" :end1 3)
(string-downcase
(concatenate 'string "u" (subseq (string type) 9)))
(string-downcase type))))
(if (and (> (length ut) 8)
(string-equal
(subseq ut (- (length ut) 9))
"long-long"))
(concatenate 'string (subseq ut 0 (- (length ut) 9)) "llong")
ut))))
*cstd-integer-types*))
;; The floating types are associated by order, so it is important that
;; order of *fp-type-mapping* and *gsl-splice-fp-types* match,
;; though the latter may be longer.
(floating-point-association *gsl-splice-fp-types*))
"Mapping the C standard types to the GSL splice name.")
(defparameter *blas-splice-fp-types*
;; Ordered by: real shortest to longest, then complex shortest to longest.
;; Must match *fp-type-mapping*.
'("s" "d" #+long-double nil
"c" "z" #+long-double nil)
"The list of floating point types that can be spliced into BLAS function names.")
(defparameter *cstd-blas-mapping*
;; The floating types are associated by order, so it is important that
;; order of *fp-type-mapping* and *blas-splice-fp-types* match,
;; though the latter may be longer.
(floating-point-association *blas-splice-fp-types*)
"Mapping the C standard types to the BLAS splice name.")
;;;;****************************************************************************
;;;; Conversions
;;;;****************************************************************************
(defun all-types (alist &optional right-side)
"A list of all types defined by symbol or definition."
(mapcar (if right-side #'rest #'first) alist))
(defun lookup-type (symbol alist &optional reverse)
"Lookup the symbol defined in the alist."
(or
(if reverse
(first (rassoc symbol alist :test #'equal))
(rest (assoc symbol alist)))
;;(error "Did not find ~a in ~a" symbol (mapcar #'first alist))
))
;;; (cl-single '(unsigned-byte 8))
;;; UNSIGNED-BYTE-8
(defun cl-single (cl-type)
"The element type name as a single symbol."
(intern (if (atom cl-type)
(princ-to-string cl-type)
(format nil "~{~a~^-~}" cl-type))
:gsl))
;;; (cl-gsl '(unsigned-byte 8))
;;; "uchar"
(defun cl-gsl (cl-type &optional prepend-underscore blas)
"The GSL splice string from the CL type."
(let ((string
(lookup-type
(lookup-type cl-type *cstd-cl-type-mapping* t)
(if blas *cstd-blas-mapping* *cstd-gsl-mapping*))))
(if (and prepend-underscore (plusp (length string)))
(concatenate 'string "_" string)
string)))
(defun cl-cffi (cl-type)
"The CFFI element type from the CL type."
(lookup-type (clean-type cl-type) *cstd-cl-type-mapping* t))
(defun cffi-cl (cffi-type)
"The CL type from the CFFI element type."
(unless (eq cffi-type :pointer)
(lookup-type cffi-type *cstd-cl-type-mapping*)))
(defun splice-name (base-name type keyword)
"Make a new C name for a data function from a base name."
(let ((pos (search keyword base-name)))
(when pos
(let ((insert (+ pos (length keyword)))
(gsltype (cl-gsl type)))
(concatenate 'string
(subseq base-name 0 insert)
(if (zerop (length gsltype)) "" "_")
(cl-gsl type)
(subseq base-name insert))))))
;;;;****************************************************************************
;;;; GSL complex types
;;;;****************************************************************************
;;; GSL defines complex numbers in a struct, and passes the struct by
;;; value. CFFI does not support call by value for structs, so we
;;; cannot use functions that call or return complex scalars.
;;; See /usr/include/gsl/gsl_complex.h
(cffi:defcstruct complex-float-c
(dat :float :count 2))
(cffi:defcstruct complex-double-c
(dat :double :count 2))
#+long-double
(cffi:defcstruct complex-long-double-c
(dat :long-double :count 2))
(defun clean-type (type)
;; SBCL at least will specify limits on the type, e.g.
;; (type-of #C(1.0 2.0))
;; (COMPLEX (DOUBLE-FLOAT 1.0 2.0))
;; This cleans that up to make
;; (clean-type (type-of #C(1.0 2.0)))
;; (COMPLEX DOUBLE-FLOAT)
(if (and (subtypep type 'complex) (listp (second type)))
(list (first type) (first (second type)))
type))
(defun complex-to-gsl (number gsl)
"Set the already-allocated GSL (foreign) struct to the CL complex number.
Returns the struct."
(cerror "Accept gibberish."
"Cannot pass complex scalars to and from GSL functions (structs passed by value).")
(let* ((cleantype (clean-type (type-of number)))
(comptype (cl-cffi (second cleantype)))
(datslot
(cffi:foreign-slot-pointer gsl (cl-cffi cleantype) 'dat)))
(setf (cffi:mem-aref datslot comptype 0) (realpart number)
(cffi:mem-aref datslot comptype 1) (imagpart number))
gsl))
(defun component-float-type (eltype)
"The float type of this type."
(if (subtypep eltype 'complex)
;; complex: use the component type
(second eltype)
eltype))
;;; Use GSL to create the complex. This actually does work on
;;; SBCL/amd64, but it shouldn't.
#+(or)
(defmfun complex-to-gsl (complex)
"gsl_complex_rect"
(((realpart complex) :double) ((imagpart complex) :double))
:c-return :pointer
:documentation
"Convert the CL complex double into a GSL complex.")
#|
(defmfun complex-abs (number)
"gsl_complex_abs"
(((complex-to-gsl number) :pointer))
:c-return :double)
(defmfun complex-arg (number)
"gsl_complex_arg"
(((complex-to-gsl number) :pointer))
:c-return :double)
|#
;;;;****************************************************************************
;;;; Common element type groups for generic functions
;;;;****************************************************************************
......
;; Conversion of numbers C->CL
;; Liam Healy, Sun May 28 2006 - 22:04
;; Time-stamp: <2008-07-06 16:33:06EDT number-conversion.lisp>
;; Time-stamp: <2009-01-13 21:40:26EST number-conversion.lisp>
;; $Id$
(in-package :gsl)
......@@ -23,7 +23,7 @@
;;;; Complex numbers
;;;;****************************************************************************
;;; GSL complex struct is defined in init/element-types.lisp.
;;; GSL complex struct is defined in init/complex-types.lisp.
(defun complex-to-cl
(gsl-complex &optional (index 0) (complex-type 'complex-double-c))
"Make a CL complex number from the GSL pointer to a complex struct or
......
;; Permissible types
;; Number types used by GSL functions, and specification conversion
;; Liam Healy 2008-12-31 21:06:34EST types.lisp
;; Time-stamp: <2009-01-08 10:28:04EST types.lisp>
;; Time-stamp: <2009-01-13 22:00:26EST types.lisp>
;; $Id: $
(in-package :gsl)
;;;;****************************************************************************
;;;; Types for CFFI (will eventually be in CFFI)
;;;; Unsigned address types size_t
;;;;****************************************************************************
(case
......@@ -19,4 +19,206 @@
(cffi:defctype sizet :uint32))
(t (error "Size of :long unrecognized")))
;;;;****************************************************************************
;;;; Type specification conversion
;;;;****************************************************************************
;; cffi-features:no-long-long doesn't work for me, but ought to be checked?
;;; The different element type forms:
;;; C standard full, or "CFFI" :unsigned-char
;;; GSL splice name "uchar"
;;; C explicit :uint8 (not used here)
;;; CL '(unsigned-byte 8)
;;; Single 'unsigned-byte-8
;;; Functions to perform conversions
;;; CL -> single in function #'cl-single
;;; CL -> GSL in function #'cl-gsl
;;; CL -> CFFI (Cstd) in function #'cl-cffi
;;; CFFI (Cstd) -> CL in function #'cffi-cl
;;; Sources of equivalence
;;; Cstd -> CL in alist *cstd-cl-type-mapping*
;;; Cstd -> GSL in alist *cstd-gsl-mapping*
;;; :long-double may be pushed onto *features* if
;;; the implementation supports long doubles in CFFI.
;;;;****************************************************************************
;;;; Basic definition
;;;;****************************************************************************
;;; Preliminary definitions. These are not used directly outside this
;;; file; they exist to define the two variables
;;; *cstd-cl-type-mapping* and *cstd-gsl-mapping* which are used by
;;; the conversion functions.
(defparameter *cstd-integer-types*
'(:char :unsigned-char
:short :unsigned-short
:int :unsigned-int
:long :unsigned-long
#+int64
:long-long
#+int64
:unsigned-long-long)
;; http://common-lisp.net/project/cffi/manual/html_node/Built_002dIn-Types.html
"List of integer types supported by CFFI, from the CFFI docs.")
(defparameter *fp-type-mapping*
'((:float . single-float) (:double . double-float)
;; For those implementations that support a separate long-double
;; type assume this mapping:
#+long-double (:long-double . long-float)
(complex-float-c . (complex single-float))
(complex-double-c . (complex double-float))
;; For those implementations that support a separate long-double
;; type assume this mapping:
#+long-double
(complex-long-double-c . (complex long-float)))
;; Ordered by: real shortest to longest, then complex shortest to longest.
"List of floating point types supported by CFFI from the CFFI docs
plus corresponding complex types.")
(defparameter *gsl-splice-int-types*
;; list | grep -i 'gsl_vector.*_alloc\b'
'("char" "int" "long" "short" "uchar" "uint" "ulong" "ushort")
"The list of integer types that can be spliced into function names.")
(defparameter *gsl-splice-fp-types*
;; list | grep -i 'gsl_vector.*_alloc\b'
;; Ordered by: real shortest to longest, then complex shortest to longest.
;; Must match *fp-type-mapping*.
'("float" "" #+long-double "long_double"
"complex_float" "complex" #+long-double "complex_long_double")
"The list of floating point types that can be spliced into function names.")
;;; Mapping alists used by conversion functions.
;;; Used by #'cl-gsl
(defparameter *cstd-cl-type-mapping*
(append
*fp-type-mapping*
(mapcar
(lambda (type)
(cons
type
(list
(if (string-equal type "uns" :end1 3)
'unsigned-byte
'signed-byte)
(* 8 (cffi:foreign-type-size type)))))
*cstd-integer-types*))
;; Be careful when reverse associating, as there may be several C
;; types that map to a single CL type.
"An alist of the C standard types as keywords, and the CL type
The exception is complex types, which don't have a definition
in the C standard; in that case, the C type is the GSL struct
definition.")
(defmacro floating-point-association (splice-list)
`(mapcar
#'cons
(mapcar #'first *fp-type-mapping*)
(subseq ,splice-list 0 (length *fp-type-mapping*))))
(defparameter *cstd-gsl-mapping*
(append
;; The integer types
(remove-if-not
(lambda (x) (find (rest x) *gsl-splice-int-types* :test 'string-equal))
(mapcar
(lambda (type)
(cons type
(let ((ut
(if (string-equal type "uns" :end1 3)
(string-downcase
(concatenate 'string "u" (subseq (string type) 9)))
(string-downcase type))))
(if (and (> (length ut) 8)
(string-equal
(subseq ut (- (length ut) 9))
"long-long"))
(concatenate 'string (subseq ut 0 (- (length ut) 9)) "llong")
ut))))
*cstd-integer-types*))
;; The floating types are associated by order, so it is important that
;; order of *fp-type-mapping* and *gsl-splice-fp-types* match,
;; though the latter may be longer.
(floating-point-association *gsl-splice-fp-types*))
"Mapping the C standard types to the GSL splice name.")
(defparameter *blas-splice-fp-types*
;; Ordered by: real shortest to longest, then complex shortest to longest.
;; Must match *fp-type-mapping*.
'("s" "d" #+long-double nil
"c" "z" #+long-double nil)
"The list of floating point types that can be spliced into BLAS function names.")
(defparameter *cstd-blas-mapping*
;; The floating types are associated by order, so it is important that
;; order of *fp-type-mapping* and *blas-splice-fp-types* match,
;; though the latter may be longer.
(floating-point-association *blas-splice-fp-types*)
"Mapping the C standard types to the BLAS splice name.")
;;;;****************************************************************************
;;;; Conversions
;;;;****************************************************************************
(defun all-types (alist &optional right-side)
"A list of all types defined by symbol or definition."
(mapcar (if right-side #'rest #'first) alist))
(defun lookup-type (symbol alist &optional reverse)
"Lookup the symbol defined in the alist."
(or
(if reverse
(first (rassoc symbol alist :test #'equal))
(rest (assoc symbol alist)))
;;(error "Did not find ~a in ~a" symbol (mapcar #'first alist))
))
;;; (cl-single '(unsigned-byte 8))
;;; UNSIGNED-BYTE-8
(defun cl-single (cl-type)
"The element type name as a single symbol."
(intern (if (atom cl-type)
(princ-to-string cl-type)
(format nil "~{~a~^-~}" cl-type))
:gsl))
;;; (cl-gsl '(unsigned-byte 8))
;;; "uchar"
(defun cl-gsl (cl-type &optional prepend-underscore blas)
"The GSL splice string from the CL type."
(let ((string
(lookup-type
(lookup-type cl-type *cstd-cl-type-mapping* t)
(if blas *cstd-blas-mapping* *cstd-gsl-mapping*))))
(if (and prepend-underscore (plusp (length string)))
(concatenate 'string "_" string)
string)))
(defun cl-cffi (cl-type)
"The CFFI element type from the CL type."
(lookup-type (clean-type cl-type) *cstd-cl-type-mapping* t))
(defun cffi-cl (cffi-type)
"The CL type from the CFFI element type."
(unless (eq cffi-type :pointer)
(lookup-type cffi-type *cstd-cl-type-mapping*)))
(defun splice-name (base-name type keyword)
"Make a new C name for a data function from a base name."
(let ((pos (search keyword base-name)))
(when pos
(let ((insert (+ pos (length keyword)))
(gsltype (cl-gsl type)))
(concatenate 'string
(subseq base-name 0 insert)
(if (zerop (length gsltype)) "" "_")
(cl-gsl type)
(subseq base-name insert))))))
......@@ -36,56 +36,56 @@
-49.08d0)))
(SCALAR -39.66d0))
(CL-ARRAY (AXPY SCALAR V1 V2)))))
(LISP-UNIT:ASSERT-ERROR 'SIMPLE-ERROR
(LET ((V1
(MAKE-MARRAY
'(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'(-34.5 8.24 3.29 -8.93
34.12 -6.15 49.27
-13.49 32.5 42.73
-17.24 43.31 -16.12
-8.25 21.44 -49.08)))
(V2
(MAKE-MARRAY
'(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'(32.5 42.73 -17.24
43.31 -16.12 -8.25
21.44 -49.08 -39.66
-49.46 19.68 -5.55
-8.82 25.37 -30.58
31.67)))
(SCALAR #C(-39.66 -49.46)))
(CL-ARRAY
(AXPY SCALAR V1 V2))))
(LISP-UNIT:ASSERT-ERROR 'SIMPLE-ERROR
(LET ((V1
(MAKE-MARRAY
'(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'(-34.5d0 8.24d0 3.29d0
-8.93d0 34.12d0
-6.15d0 49.27d0
-13.49d0 32.5d0
42.73d0 -17.24d0
43.31d0 -16.12d0
-8.25d0 21.44d0
-49.08d0)))
(V2
(MAKE-MARRAY
'(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'(32.5d0 42.73d0
-17.24d0 43.31d0
-16.12d0 -8.25d0
21.44d0 -49.08d0
-39.66d0 -49.46d0
19.68d0 -5.55d0
-8.82d0 25.37d0
-30.58d0 31.67d0)))
(SCALAR
#C(-39.66d0 -49.46d0)))
(CL-ARRAY
(AXPY SCALAR V1 V2)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST
#(#C(1808.3204 1422.3015) #C(-589.3992 234.75043)
#C(-1673.498 -1451.916) #C(-2599.8237 -1950.9608)
#C(784.81586 -3351.5815) #C(2845.531 -870.5343)
#C(222.45422 1149.8602) #C(-3308.3872 917.76044)))
(MULTIPLE-VALUE-LIST
(LET ((V1
(MAKE-MARRAY '(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'(-34.5 8.24 3.29 -8.93 34.12
-6.15 49.27 -13.49 32.5 42.73
-17.24 43.31 -16.12 -8.25 21.44
-49.08)))
(V2
(MAKE-MARRAY '(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'(32.5 42.73 -17.24 43.31 -16.12
-8.25 21.44 -49.08 -39.66 -49.46
19.68 -5.55 -8.82 25.37 -30.58
31.67)))
(SCALAR #C(-39.66 -49.46)))
(CL-ARRAY (AXPY SCALAR V1 V2)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST
#(#C(1808.3204d0 1422.3016000000002d0)
#C(-589.3992d0 234.75039999999998d0)
#C(-1673.4981999999998d0 -1451.9162000000001d0)
#C(-2599.8236d0 -1950.9608000000003d0)
#C(784.8158000000002d0 -3351.5818d0)
#C(2845.5309999999995d0 -870.5342d0)
#C(222.45420000000001d0 1149.8601999999998d0)
#C(-3308.3871999999997d0 917.7603999999995d0)))
(MULTIPLE-VALUE-LIST
(LET ((V1
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'(-34.5d0 8.24d0 3.29d0 -8.93d0
34.12d0 -6.15d0 49.27d0 -13.49d0
32.5d0 42.73d0 -17.24d0 43.31d0
-16.12d0 -8.25d0 21.44d0
-49.08d0)))
(V2
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'(32.5d0 42.73d0 -17.24d0 43.31d0
-16.12d0 -8.25d0 21.44d0
-49.08d0 -39.66d0 -49.46d0
19.68d0 -5.55d0 -8.82d0 25.37d0
-30.58d0 31.67d0)))
(SCALAR #C(-39.66d0 -49.46d0)))
(CL-ARRAY (AXPY SCALAR V1 V2))))))
......@@ -42,63 +42,67 @@
(-49.08d0 -39.66d0 -49.46d0))))
(S1 19.68d0))
(CL-ARRAY (INVERSE-MATRIX-PRODUCT M1 M2 S1)))))
(LISP-UNIT:ASSERT-ERROR 'SIMPLE-ERROR
(LET ((M1
(MAKE-MARRAY
'(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5 8.24 3.29 -8.93
34.12 -6.15)
(-8.93 34.12 -6.15
49.27 -13.49 32.5)
(49.27 -13.49 32.5
42.73 -17.24
43.31))))
(M2
(MAKE-MARRAY
'(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'((42.73 -17.24 43.31
-16.12 -8.25 21.44)
(-16.12 -8.25 21.44
-49.08 -39.66 -49.46)
(-49.08 -39.66 -49.46
19.68 -5.55 -8.82))))
(S1 #C(19.68 -5.55)))
(CL-ARRAY
(INVERSE-MATRIX-PRODUCT M1 M2
S1))))
(LISP-UNIT:ASSERT-ERROR 'SIMPLE-ERROR
(LET ((M1
(MAKE-MARRAY
'(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5d0 8.24d0 3.29d0
-8.93d0 34.12d0
-6.15d0)
(-8.93d0 34.12d0
-6.15d0 49.27d0
-13.49d0 32.5d0)
(49.27d0 -13.49d0
32.5d0 42.73d0
-17.24d0 43.31d0))))
(M2
(MAKE-MARRAY
'(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((42.73d0 -17.24d0
43.31d0 -16.12d0
-8.25d0 21.44d0)
(-16.12d0 -8.25d0
21.44d0 -49.08d0
-39.66d0 -49.46d0)
(-49.08d0 -39.66d0
-49.46d0 19.68d0
-5.55d0 -8.82d0))))
(S1 #C(19.68d0 -5.55d0)))
(CL-ARRAY
(INVERSE-MATRIX-PRODUCT M1 M2
S1))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST
#2A((#C(-28.191956 35.052914) #C(-13.378113 28.736774)
#C(5.4411864 -3.0717967))
(#C(5.1555862 -11.330741) #C(-33.49236 -12.241006)
#C(-10.527424 20.3357))
(#C(-0.71776515 27.669863) #C(20.046593 11.972884)
#C(-1.5907472 4.2853727))))
(MULTIPLE-VALUE-LIST
(LET ((M1
(MAKE-MARRAY '(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5 8.24 3.29 -8.93 34.12
-6.15)
(-8.93 34.12 -6.15 49.27 -13.49
32.5)
(49.27 -13.49 32.5 42.73 -17.24
43.31))))
(M2
(MAKE-MARRAY '(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'((42.73 -17.24 43.31 -16.12 -8.25
21.44)
(-16.12 -8.25 21.44 -49.08
-39.66 -49.46)
(-49.08 -39.66 -49.46 19.68
-5.55 -8.82))))
(S1 #C(19.68 -5.55)))
(CL-ARRAY (INVERSE-MATRIX-PRODUCT M1 M2 S1)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST
#2A((#C(-28.191956546692666d0 35.052909000810466d0)
#C(-13.378112744466048d0 28.736769449657892d0)
#C(5.441186078394205d0 -3.071795512940909d0))
(#C(5.155586147762611d0 -11.330740389402608d0)
#C(-33.49236130612107d0 -12.241006594989193d0)
#C(-10.527425621635187d0 20.33570209413866d0))
(#C(-0.7177649743298794d0 27.669860728641122d0)
#C(20.046591997869093d0 11.972882797431007d0)
#C(-1.590747546093172d0 4.285372608973592d0))))
(MULTIPLE-VALUE-LIST
(LET ((M1
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5d0 8.24d0 3.29d0 -8.93d0
34.12d0 -6.15d0)
(-8.93d0 34.12d0 -6.15d0 49.27d0
-13.49d0 32.5d0)
(49.27d0 -13.49d0 32.5d0 42.73d0
-17.24d0 43.31d0))))
(M2
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((42.73d0 -17.24d0 43.31d0
-16.12d0 -8.25d0 21.44d0)
(-16.12d0 -8.25d0 21.44d0
-49.08d0 -39.66d0 -49.46d0)
(-49.08d0 -39.66d0 -49.46d0
19.68d0 -5.55d0 -8.82d0))))
(S1 #C(19.68d0 -5.55d0)))
(CL-ARRAY (INVERSE-MATRIX-PRODUCT M1 M2 S1)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(-1.1747805 -0.26507667 1.3326154))
(MULTIPLE-VALUE-LIST
......
This diff is collapsed.
......@@ -42,65 +42,68 @@
(-49.08d0 -39.66d0 -49.46d0))))
(S1 19.68d0))
(CL-ARRAY (MATRIX-PRODUCT-TRIANGULAR M1 M2 S1)))))
(LISP-UNIT:ASSERT-ERROR 'SIMPLE-ERROR
(LET ((M1
(MAKE-MARRAY
'(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5 8.24 3.29 -8.93
34.12 -6.15)
(-8.93 34.12 -6.15
49.27 -13.49 32.5)
(49.27 -13.49 32.5
42.73 -17.24
43.31))))
(M2
(MAKE-MARRAY
'(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'((42.73 -17.24 43.31
-16.12 -8.25 21.44)
(-16.12 -8.25 21.44
-49.08 -39.66 -49.46)
(-49.08 -39.66 -49.46
19.68 -5.55 -8.82))))
(S1 #C(19.68 -5.55)))
(CL-ARRAY
(MATRIX-PRODUCT-TRIANGULAR M1
M2
S1))))
(LISP-UNIT:ASSERT-ERROR 'SIMPLE-ERROR
(LET ((M1
(MAKE-MARRAY
'(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5d0 8.24d0 3.29d0
-8.93d0 34.12d0
-6.15d0)
(-8.93d0 34.12d0
-6.15d0 49.27d0
-13.49d0 32.5d0)
(49.27d0 -13.49d0
32.5d0 42.73d0
-17.24d0 43.31d0))))
(M2
(MAKE-MARRAY
'(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((42.73d0 -17.24d0
43.31d0 -16.12d0
-8.25d0 21.44d0)
(-16.12d0 -8.25d0
21.44d0 -49.08d0
-39.66d0 -49.46d0)
(-49.08d0 -39.66d0
-49.46d0 19.68d0
-5.55d0 -8.82d0))))
(S1 #C(19.68d0 -5.55d0)))
(CL-ARRAY
(MATRIX-PRODUCT-TRIANGULAR M1
M2
S1))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST
#2A((#C(-66397.9 18986.908) #C(-56335.145 48514.305)
#C(-18830.469 -13449.601))
(#C(38337.09 -49128.918) #C(42681.344 -22972.445)
#C(50375.406 -50562.54))
(#C(42453.223 -42606.086)
#C(-13764.868 -48835.813)
#C(8910.526 -4389.117))))
(MULTIPLE-VALUE-LIST
(LET ((M1
(MAKE-MARRAY '(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5 8.24 3.29 -8.93 34.12
-6.15)
(-8.93 34.12 -6.15 49.27 -13.49
32.5)
(49.27 -13.49 32.5 42.73 -17.24
43.31))))
(M2
(MAKE-MARRAY '(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'((42.73 -17.24 43.31 -16.12 -8.25
21.44)
(-16.12 -8.25 21.44 -49.08
-39.66 -49.46)
(-49.08 -39.66 -49.46 19.68
-5.55 -8.82))))
(S1 #C(19.68 -5.55)))
(CL-ARRAY (MATRIX-PRODUCT-TRIANGULAR M1 M2 S1)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST
#2A((#C(-66397.89753899998d0 18986.908142999997d0)
#C(-56335.14258600001d0 48514.306278000004d0)
#C(-18830.468709d0 -13449.603000000003d0))
(#C(38337.08717099999d0 -49128.917505000005d0)
#C(42681.34176d0 -22972.447481999992d0)
#C(50375.404416000005d0 -50562.535017d0))
(#C(42453.21956399999d0 -42606.08134200001d0)
#C(-13764.866562000001d0 -48835.809623999994d0)
#C(8910.526581d0 -4389.116526d0))))
(MULTIPLE-VALUE-LIST
(LET ((M1
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5d0 8.24d0 3.29d0 -8.93d0
34.12d0 -6.15d0)
(-8.93d0 34.12d0 -6.15d0 49.27d0
-13.49d0 32.5d0)
(49.27d0 -13.49d0 32.5d0 42.73d0
-17.24d0 43.31d0))))
(M2
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((42.73d0 -17.24d0 43.31d0
-16.12d0 -8.25d0 21.44d0)
(-16.12d0 -8.25d0 21.44d0
-49.08d0 -39.66d0 -49.46d0)
(-49.08d0 -39.66d0 -49.46d0
19.68d0 -5.55d0 -8.82d0))))
(S1 #C(19.68d0 -5.55d0)))
(CL-ARRAY (MATRIX-PRODUCT-TRIANGULAR M1 M2 S1)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(-1473.7527 -854.58527 1407.5751))
(MULTIPLE-VALUE-LIST
......
......@@ -54,89 +54,87 @@
(S1 -41.67d0)
(S2 42.0d0))
(CL-ARRAY (MATRIX-PRODUCT M1 M2 M3 S1 S2)))))
(LISP-UNIT:ASSERT-ERROR 'SIMPLE-ERROR
(LET ((M1
(MAKE-MARRAY
'(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5 8.24 3.29 -8.93
34.12 -6.15)
(-8.93 34.12 -6.15
49.27 -13.49 32.5)
(49.27 -13.49 32.5
42.73 -17.24
43.31))))
(M2
(MAKE-MARRAY
'(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'((42.73 -17.24 43.31
-16.12 -8.25 21.44)
(-16.12 -8.25 21.44
-49.08 -39.66 -49.46)
(-49.08 -39.66 -49.46
19.68 -5.55 -8.82))))
(M3
(MAKE-MARRAY
'(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'((19.68 -5.55 -8.82
25.37 -30.58 31.67)
(25.37 -30.58 31.67
29.36 -33.24 -27.03)
(29.36 -33.24 -27.03
-41.67 42.0
-20.81))))
(S1 #C(-41.67 42.0))
(S2 #C(42.0 -20.81)))
(CL-ARRAY
(MATRIX-PRODUCT M1 M2 M3 S1
S2))))
(LISP-UNIT:ASSERT-ERROR 'SIMPLE-ERROR
(LET ((M1
(MAKE-MARRAY
'(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5d0 8.24d0 3.29d0
-8.93d0 34.12d0
-6.15d0)
(-8.93d0 34.12d0
-6.15d0 49.27d0
-13.49d0 32.5d0)
(49.27d0 -13.49d0
32.5d0 42.73d0
-17.24d0 43.31d0))))
(M2
(MAKE-MARRAY
'(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((42.73d0 -17.24d0
43.31d0 -16.12d0
-8.25d0 21.44d0)
(-16.12d0 -8.25d0
21.44d0 -49.08d0
-39.66d0 -49.46d0)
(-49.08d0 -39.66d0
-49.46d0 19.68d0
-5.55d0 -8.82d0))))
(M3
(MAKE-MARRAY
'(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((19.68d0 -5.55d0
-8.82d0 25.37d0
-30.58d0 31.67d0)
(25.37d0 -30.58d0
31.67d0 29.36d0
-33.24d0 -27.03d0)
(29.36d0 -33.24d0
-27.03d0 -41.67d0
42.0d0 -20.81d0))))
(S1 #C(-41.67d0 42.0d0))
(S2 #C(42.0d0 -20.81d0)))
(CL-ARRAY
(MATRIX-PRODUCT M1 M2 M3 S1
S2))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST
#2A((#C(140927.64 -143005.16) #C(72986.305 -201162.86)
#C(65956.16 9034.215))
(#C(-102500.125 118033.82) #C(-147775.0 58489.473)
#C(-9582.236 190719.98))
(#C(-16551.412 336505.22) #C(-7668.242 385798.25)
#C(45313.527 140641.6))))
(MULTIPLE-VALUE-LIST
(LET ((M1
(MAKE-MARRAY '(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5 8.24 3.29 -8.93 34.12
-6.15)
(-8.93 34.12 -6.15 49.27 -13.49
32.5)
(49.27 -13.49 32.5 42.73 -17.24
43.31))))
(M2
(MAKE-MARRAY '(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'((42.73 -17.24 43.31 -16.12 -8.25
21.44)
(-16.12 -8.25 21.44 -49.08
-39.66 -49.46)
(-49.08 -39.66 -49.46 19.68
-5.55 -8.82))))
(M3
(MAKE-MARRAY '(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'((19.68 -5.55 -8.82 25.37 -30.58
31.67)
(25.37 -30.58 31.67 29.36 -33.24
-27.03)
(29.36 -33.24 -27.03 -41.67 42.0
-20.81))))
(S1 #C(-41.67 42.0))
(S2 #C(42.0 -20.81)))
(CL-ARRAY (MATRIX-PRODUCT M1 M2 M3 S1 S2)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST
#2A((#C(140927.638311d0 -143005.14965699997d0)
#C(72986.314354d0 -201162.870342d0)
#C(65956.156396d0 9034.212785d0))
(#C(-102500.143082d0 118033.81826900001d0)
#C(-147775.01365700003d0 58489.469011999994d0)
#C(-9582.23530800004d0 190719.98257199998d0))
(#C(-16551.43698600002d0 336505.21680500003d0)
#C(-7668.248943000013d0 385798.22976300004d0)
#C(45313.51313299998d0 140641.59851399998d0))))
(MULTIPLE-VALUE-LIST
(LET ((M1
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5d0 8.24d0 3.29d0 -8.93d0
34.12d0 -6.15d0)
(-8.93d0 34.12d0 -6.15d0 49.27d0
-13.49d0 32.5d0)
(49.27d0 -13.49d0 32.5d0 42.73d0
-17.24d0 43.31d0))))
(M2
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((42.73d0 -17.24d0 43.31d0
-16.12d0 -8.25d0 21.44d0)
(-16.12d0 -8.25d0 21.44d0
-49.08d0 -39.66d0 -49.46d0)
(-49.08d0 -39.66d0 -49.46d0
19.68d0 -5.55d0 -8.82d0))))
(M3
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((19.68d0 -5.55d0 -8.82d0 25.37d0
-30.58d0 31.67d0)
(25.37d0 -30.58d0 31.67d0
29.36d0 -33.24d0 -27.03d0)
(29.36d0 -33.24d0 -27.03d0
-41.67d0 42.0d0 -20.81d0))))
(S1 #C(-41.67d0 42.0d0))
(S2 #C(42.0d0 -20.81d0)))
(CL-ARRAY (MATRIX-PRODUCT M1 M2 M3 S1 S2)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(72971.1 60998.137 -184676.98))
(MULTIPLE-VALUE-LIST
......@@ -173,67 +171,59 @@
(S1 -49.08d0)
(S2 -39.66d0))
(CL-ARRAY (MATRIX-PRODUCT M1 V1 V2 S1 S2)))))
(LISP-UNIT:ASSERT-ERROR 'SIMPLE-ERROR
(LET ((M1
(MAKE-MARRAY
'(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5 8.24 3.29 -8.93
34.12 -6.15)
(-8.93 34.12 -6.15
49.27 -13.49 32.5)
(49.27 -13.49 32.5
42.73 -17.24
43.31))))
(V1
(MAKE-MARRAY
'(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'(42.73 -17.24 43.31
-16.12 -8.25 21.44)))
(V2
(MAKE-MARRAY
'(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'(-16.12 -8.25 21.44
-49.08 -39.66
-49.46)))
(S1 #C(-49.08 -39.66))
(S2 #C(-39.66 -49.46)))
(CL-ARRAY
(MATRIX-PRODUCT M1 V1 V2 S1
S2))))
(LISP-UNIT:ASSERT-ERROR 'SIMPLE-ERROR
(LET ((M1
(MAKE-MARRAY
'(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5d0 8.24d0 3.29d0
-8.93d0 34.12d0
-6.15d0)
(-8.93d0 34.12d0
-6.15d0 49.27d0
-13.49d0 32.5d0)
(49.27d0 -13.49d0
32.5d0 42.73d0
-17.24d0 43.31d0))))
(V1
(MAKE-MARRAY
'(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'(42.73d0 -17.24d0
43.31d0 -16.12d0
-8.25d0 21.44d0)))
(V2
(MAKE-MARRAY
'(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'(-16.12d0 -8.25d0
21.44d0 -49.08d0
-39.66d0 -49.46d0)))
(S1 #C(-49.08d0 -39.66d0))
(S2 #C(-39.66d0 -49.46d0)))
(CL-ARRAY
(MATRIX-PRODUCT M1 V1 V2 S1
S2)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST
#(#C(124163.57 -3332.2988) #C(119793.48 -166378.0)
#C(-189845.23 -81764.12)))
(MULTIPLE-VALUE-LIST
(LET ((M1
(MAKE-MARRAY '(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5 8.24 3.29 -8.93 34.12
-6.15)
(-8.93 34.12 -6.15 49.27 -13.49
32.5)
(49.27 -13.49 32.5 42.73 -17.24
43.31))))
(V1
(MAKE-MARRAY '(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'(42.73 -17.24 43.31 -16.12 -8.25
21.44)))
(V2
(MAKE-MARRAY '(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'(-16.12 -8.25 21.44 -49.08 -39.66
-49.46)))
(S1 #C(-49.08 -39.66))
(S2 #C(-39.66 -49.46)))
(CL-ARRAY (MATRIX-PRODUCT M1 V1 V2 S1 S2)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST
#(#C(124163.58013199999d0 -3332.2974459999896d0)
#C(119793.47618999999d0 -166378.00423d0)
#C(-189845.21679399998d0 -81764.10481799999d0)))
(MULTIPLE-VALUE-LIST
(LET ((M1
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5d0 8.24d0 3.29d0 -8.93d0
34.12d0 -6.15d0)
(-8.93d0 34.12d0 -6.15d0 49.27d0
-13.49d0 32.5d0)
(49.27d0 -13.49d0 32.5d0 42.73d0
-17.24d0 43.31d0))))
(V1
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'(42.73d0 -17.24d0 43.31d0
-16.12d0 -8.25d0 21.44d0)))
(V2
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'(-16.12d0 -8.25d0 21.44d0
-49.08d0 -39.66d0 -49.46d0)))
(S1 #C(-49.08d0 -39.66d0))
(S2 #C(-39.66d0 -49.46d0)))
(CL-ARRAY (MATRIX-PRODUCT M1 V1 V2 S1 S2))))))
......@@ -62,32 +62,41 @@
-13.49d0)))
(SCALAR 32.5d0))
(CL-ARRAY (SCALE SCALAR V1)))))
(LISP-UNIT:ASSERT-ERROR 'SIMPLE-ERROR
(LET ((V1
(MAKE-MARRAY
'(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'(-34.5 8.24 3.29 -8.93
34.12 -6.15 49.27
-13.49 32.5 42.73
-17.24 43.31 -16.12
-8.25 21.44 -49.08)))
(SCALAR #C(32.5 42.73)))
(CL-ARRAY (SCALE SCALAR V1))))
(LISP-UNIT:ASSERT-ERROR 'SIMPLE-ERROR
(LET ((V1
(MAKE-MARRAY
'(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'(-34.5d0 8.24d0 3.29d0
-8.93d0 34.12d0
-6.15d0 49.27d0
-13.49d0 32.5d0
42.73d0 -17.24d0
43.31d0 -16.12d0
-8.25d0 21.44d0
-49.08d0)))
(SCALAR
#C(32.5d0 42.73d0)))
(CL-ARRAY (SCALE SCALAR V1)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST
#(#C(-1473.3452 -1206.385) #C(488.5039 -149.64331)
#C(1371.6895 1258.0725) #C(2177.7026 1666.8821)
#C(-769.6029 2777.45) #C(-2410.9363 670.9099)
#C(-171.37753 -956.9326) #C(2793.9885 -678.9689)))
(MULTIPLE-VALUE-LIST
(LET ((V1
(MAKE-MARRAY '(COMPLEX SINGLE-FLOAT)
:INITIAL-CONTENTS
'(-34.5 8.24 3.29 -8.93 34.12
-6.15 49.27 -13.49 32.5 42.73
-17.24 43.31 -16.12 -8.25 21.44
-49.08)))
(SCALAR #C(32.5 42.73)))
(CL-ARRAY (SCALE SCALAR V1)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST
#(#C(-1473.3452d0 -1206.385d0)
#C(488.5039d0 -149.64329999999998d0)
#C(1371.6895d0 1258.0725999999997d0)
#C(2177.7027d0 1666.8821d0)
#C(-769.6028999999996d0 2777.45d0)
#C(-2410.9363d0 670.9098000000001d0)
#C(-171.3775d0 -956.9326d0)
#C(2793.9883999999997d0 -678.9687999999999d0)))
(MULTIPLE-VALUE-LIST
(LET ((V1
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'(-34.5d0 8.24d0 3.29d0 -8.93d0
34.12d0 -6.15d0 49.27d0 -13.49d0
32.5d0 42.73d0 -17.24d0 43.31d0
-16.12d0 -8.25d0 21.44d0
-49.08d0)))
(SCALAR #C(32.5d0 42.73d0)))
(CL-ARRAY (SCALE SCALAR V1))))))
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