Newer
Older
;; BLAS level 2, Matrix-vector operations
;; Liam Healy, Wed Apr 26 2006 - 21:08
;; Time-stamp: <2011-02-19 13:08:23EST blas2.lisp>
;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program. If not, see <http://www.gnu.org/licenses/>.
;;; /usr/include/gsl/gsl_blas.h
;;;;****************************************************************************
;;;; Options
;;;;****************************************************************************
(fsbv:defcenum-aux cblas-transpose)
(fsbv:defcenum-aux cblas-uplo)
(fsbv:defcenum-aux cblas-diag)
;;;;****************************************************************************
;;;;****************************************************************************
(defun matrix-product-dimensions (a b)
(if (typep b 'grid:matrix)
(list (first (grid:dimensions a))
(second (grid:dimensions b)))
(first (grid:dimensions a))))
((A grid:matrix) (x vector)
y
(alpha 1) (beta 1) (TransA :notrans) TransB
&aux
(grid:ensure-foreign-array
y (matrix-product-dimensions A x) nil element-type 0)))
("gsl_blas_" :type "gemv")
((transa cblas-transpose) (alpha :element-c-type) ((mpointer A) :pointer)
((mpointer x) :pointer) (beta :element-c-type) ((mpointer yarr) :pointer))
"If the second and third arguments are vectors, compute
the matrix-vector product and sum
for TransA = :notrans, :trans, :conjtrans.
If the second and third arguments are matrices, compute
the matrix-matrix product and sum C = alpha
op(A) op(B) + beta C where op(A) = A, A^T, A^H for TransA =
:notrans, :trans, :conjtrans and similarly for the
parameter TransB.")
((A grid:matrix) (x vector)
&optional (alpha 1) (uplo :upper) (TransA :notrans)
(diag :nonunit) (side :left))
("gsl_blas_" :type "trmv")
((uplo cblas-uplo) (TransA cblas-transpose) (diag cblas-diag)
((mpointer A) :pointer) ((mpointer x) :pointer))
:definition :generic
:element-types :float-complex
"If the second argument is a vector, compute
for the triangular matrix A, where op(A) = A, A^T, A^H for
TransA = :NoTrans, :Trans, :ConjTrans. When Uplo
is :Upper then the upper triangle of A is used, and when
Uplo is :Lower then the lower triangle of A is used. If
Diag is :NonUnit then the diagonal of the matrix is used,
but if Diag is :Unit then the diagonal elements of the
matrix A are taken as unity and are not referenced.
If the second argument is a matrix, compute
the matrix-matrix product B = alpha op(A) B
:Right. The matrix A is triangular and op(A) = A, A^T, A^H
for TransA = :NoTrans, :Trans, :ConjTrans When Uplo
is :Upper then the upper triangle of A is used, and when
Uplo is :Lower then the lower triangle of A is used. If
Diag is :NonUnit then the diagonal of A is used, but if
Diag is :Unit then the diagonal elements of the matrix A
are taken as unity and are not referenced.")
((A grid:matrix) (x vector)
&optional (alpha 1) (uplo :upper) (TransA :notrans)
(diag :nonunit) (side :left))
("gsl_blas_" :type "trsv")
((uplo cblas-uplo) (TransA cblas-transpose) (diag cblas-diag)
((mpointer A) :pointer) ((mpointer x) :pointer))
:definition :generic
:element-types :float-complex
:inputs (A x)
:outputs (x)
:documentation ; FDL
"If the second argument is a vector, compute
inv(op(A)) x for x, where op(A) = A, A^T, A^H for
TransA = :NoTrans, :Trans, :ConjTrans. When Uplo is
:Upper then the upper triangle of A is used, and when Uplo is
:Lower then the lower triangle of A is used. If Diag is
:NonUnit then the diagonal of the matrix is used, but if Diag
is :Unit then the diagonal elements of the matrix A are taken
as unity and are not referenced.
If the second argument is a matrix, compute
the inverse-matrix matrix product B = alpha op(inv(A))B if
Side is :Left and B = alpha B op(inv(A)) if Side is
:Right. The matrix A is triangular and op(A) = A, A^T, A^H
for TransA = :NoTrans, :Trans, :ConjTrans When
Uplo is :Upper then the upper triangle of A is used, and
when Uplo is :Lower then the lower triangle of A is
used. If Diag is :NonUnit then the diagonal of A is used,
but if Diag is :Unit then the diagonal elements of the
matrix A are taken as unity and are not referenced.")
(defmfun matrix-product-symmetric
((A grid:matrix) (x vector)
&optional y
(alpha 1) (beta 1) (uplo :upper) (side :left)
&aux
(yarr
(or y
(grid:make-foreign-array element-type :dimensions (matrix-product-dimensions A x)
:initial-element 0))))
("gsl_blas_" :type "symv")
((uplo cblas-uplo) (alpha :element-c-type) ((mpointer A) :pointer)
((mpointer x) :pointer) (beta :element-c-type) ((mpointer yarr) :pointer))
:inputs (A x yarr)
:outputs (yarr)
"If the second and third arguments are vectors, compute
the matrix-vector product and sum y = alpha A
x + beta y for the symmetric matrix A. Since the matrix A is
symmetric only its upper half or lower half need to be
stored. When Uplo is :Upper then the upper triangle and
diagonal of A are used, and when Uplo is :Lower then the
lower triangle and diagonal of A are used.
If the second and third arguments are matrices, compute
the matrix-matrix product and sum C = alpha A
B + beta C for Side is :Left and C = alpha B A + beta C
for Side is :Right, where the matrix A is symmetric. When
Uplo is :Upper then the upper triangle and diagonal of A
are used, and when Uplo is :Lower then the lower triangle
and diagonal of A are used.")
((A grid:matrix) (x vector)
(y (grid:make-foreign-array element-type :dimensions (matrix-product-dimensions A x)
(alpha 1) (beta 1) (uplo :upper) (side :left))
("gsl_blas_" :type "hemv")
((uplo cblas-uplo) (alpha :element-c-type) ((mpointer A) :pointer)
((mpointer x) :pointer) (beta :element-c-type) ((mpointer y) :pointer))
:definition :generic
:element-types :complex
:inputs (A x y)
:outputs (y)
:documentation ; FDL
"If the second and third arguments are vectors, compute the
matrix-vector product and sum y = alpha A x + beta y for the
hermitian matrix A. Since the matrix A is hermitian only its upper
half or lower half need to be stored. When Uplo is :upper then
the upper triangle and diagonal of A are used, and when Uplo is
:lower then the lower triangle and diagonal of A are used. The
imaginary elements of the diagonal are automatically assumed to be
zero and are not referenced. If the second and third arguments are
matrices, compute the matrix-matrix product and sum C = alpha A B +
beta C if Side is :left and C = \alpha B A + \beta C if Side
is :right, where the matrix A is hermitian. When Uplo is
:upper then the upper triangle and diagonal of A are used, and
when Uplo is :lower then the lower triangle and diagonal of A
are used. The imaginary elements of the diagonal are automatically
set to zero.")
(defmfun rank-1-update (alpha (x vector) (y vector) (A grid:matrix))
((alpha :element-c-type) ((mpointer x) :pointer)
((mpointer y) :pointer) ((mpointer A) :pointer))
:inputs (x y A)
:outputs (A)
:documentation ; FDL
"The rank-1 update A = alpha x y^T + A of the matrix A.")
(defmfun conjugate-rank-1-update (alpha (x vector) (y vector) (A grid:matrix))
((alpha :element-c-type) ((mpointer x) :pointer)
((mpointer y) :pointer) ((mpointer A) :pointer))
:inputs (x y A)
:outputs (A)
:documentation ; FDL
"The conjugate rank-1 update A = alpha x y^H + A of the matrix A.")
((x vector) (A grid:matrix)
&optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
("gsl_blas_" :type "syr")
((uplo cblas-uplo) (alpha :element-c-type)
((mpointer x) :pointer) ((mpointer A) :pointer))
:definition :generic
:element-types :float
:inputs (x A)
:outputs (A)
:documentation ; FDL
"If the first argument is a vector,
the symmetric rank-1 update A = \alpha x x^T + A of the symmetric
matrix A. Since the matrix A is symmetric only its upper half or
lower half need to be stored. When Uplo is :Upper then the upper
triangle and diagonal of A are used, and when Uplo is :Lower then
the lower triangle and diagonal of A are used. If the first
argument is a matrix, a rank-k update of the symmetric matrix C, C =
\alpha A A^T + \beta C when Trans is CblasNoTrans and C = \alpha A^T
A + \beta C when Trans is CblasTrans. Since the matrix C is
symmetric only its upper half or lower half need to be stored. When
Uplo is CblasUpper then the upper triangle and diagonal of C are
used, and when Uplo is CblasLower then the lower triangle and
diagonal of C are used.")
((x vector) (A grid:matrix)
&optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
("gsl_blas_" :type "her")
((uplo cblas-uplo) (alpha :element-c-type)
((mpointer x) :pointer) ((mpointer A) :pointer))
:definition :generic
:element-types :complex
:inputs (x A)
:outputs (A)
:documentation ; FDL
"If the first argument is a vector,
compute the hermitian rank-1 update A = alpha x x^H + A of the
hermitian matrix A. Since the matrix A is hermitian only its upper
half or lower half need to be stored. When Uplo is :upper then the
upper triangle and diagonal of A are used, and when Uplo is
:lower then the lower triangle and diagonal of A are used. The
imaginary elements of the diagonal are automatically set to zero.
If the first argument is a matrix, compute a rank-k update of the
hermitian matrix C, C = \alpha A A^H + \beta C when Trans is
:notrans and C = \alpha A^H A + \beta C when Trans is
:trans. Since the matrix C is hermitian only its upper half or
lower half need to be stored. When Uplo is :upper then the upper
triangle and diagonal of C are used, and when Uplo is :lower
then the lower triangle and diagonal of C are used. The imaginary
elements of the diagonal are automatically set to zero.")
((x vector) (y vector) (A grid:matrix)
&optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
("gsl_blas_" :type "syr2")
((uplo cblas-uplo) (alpha :element-c-type)
((mpointer x) :pointer) ((mpointer y) :pointer)
((mpointer A) :pointer))
:definition :generic
:element-types :float
:inputs (x y A)
:outputs (A)
:documentation ; FDL
"If the first two arguments are vectors, compute
the symmetric rank-2 update A = alpha x y^T + alpha y x^T + A of the
symmetric matrix A. Since the matrix A is symmetric only its upper
half or lower half need to be stored. When Uplo is :upper then the
upper triangle and diagonal of A are used, and when Uplo is :lower
then the lower triangle and diagonal of A are used. If the first
two arguments are matrices, compute a rank-2k update of the
symmetric matrix C, C = \alpha A B^T + \alpha B A^T + \beta C when
Trans is :notrans and C = \alpha A^T B + \alpha B^T A + \beta C
when Trans is :trans. Since the matrix C is symmetric only its
upper half or lower half need to be stored. When Uplo is :upper
then the upper triangle and diagonal of C are used, and when Uplo is
:lower then the lower triangle and diagonal of C are used.")
((x vector) (y vector) (A grid:matrix)
&optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
("gsl_blas_" :type "her2")
((uplo cblas-uplo) (alpha :element-c-type)
((mpointer x) :pointer) ((mpointer A) :pointer))
:definition :generic
:element-types :complex
:inputs (x A)
:outputs (A)
:documentation ; FDL
"If the first two arguments are vectors, compute the
hermitian rank-2 update A = alpha x y^H + alpha^* y x^H A of
the hermitian matrix A. Since the matrix A is hermitian only its
upper half or lower half need to be stored. When uplo is :upper
then the upper triangle and diagonal of A are used, and when uplo is
:lower then the lower triangle and diagonal of A are used. The
imaginary elements of the diagonal are automatically set to zero.
If the first two arguments are matrices, compute a rank-2k update of
the hermitian matrix C, C = \alpha A B^H + \alpha^* B A^H + \beta C
when Trans is :notrans and C = \alpha A^H B + \alpha^* B^H A +
\beta C when Trans is :conjtrans. Since the matrix C is
hermitian only its upper half or lower half need to be stored. When
Uplo is :upper then the upper triangle and diagonal of C are
used, and when Uplo is :lower then the lower triangle and
diagonal of C are used. The imaginary elements of the diagonal are
automatically set to zero.")
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
(generate-all-array-tests matrix-product #+fsbv :float-complex #-fsbv :float
(let ((m1 (array-default '(3 3)))
(v1 (array-default 3))
(v2 (array-default 3))
(s1 (scalar-default))
(s2 (scalar-default)))
(grid:copy-to (matrix-product m1 v1 v2 s1 s2))))
(generate-all-array-tests matrix-product-triangular
#+fsbv :float-complex #-fsbv :float
(let ((m1 (array-default '(3 3)))
(v1 (array-default 3))
(s1 (scalar-default)))
(grid:copy-to (matrix-product-triangular m1 v1 s1))))
(generate-all-array-tests inverse-matrix-product
#+fsbv :float-complex #-fsbv :float
(let ((m1 (array-default '(3 3)))
(v1 (array-default 3))
(s1 (scalar-default)))
(grid:copy-to (inverse-matrix-product m1 v1 s1))))
(generate-all-array-tests matrix-product-symmetric :float
(let ((m1 (array-default '(3 3)))
(v3 (array-default 3))
(grid:copy-to (matrix-product-symmetric m1 v1 v3 s1 s2))))
(generate-all-array-tests matrix-product-hermitian :complex
(let ((m1 (array-default '(3 3)))
(v1 (array-default 3))
(v2 (array-default 3))
(s1 (scalar-default))
(s2 (scalar-default)))
(grid:copy-to (matrix-product-hermitian m1 v1 v2 s1 s2))))
(generate-all-array-tests rank-1-update :float-complex
(let ((m1 (array-default '(3 3)))
(v1 (array-default 3))
(v2 (array-default 3))
(s1 (scalar-default)))
(grid:copy-to (rank-1-update s1 v1 v2 m1))))