;; 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/>. ;; $Id$ (in-package :gsl) ;;; /usr/include/gsl/gsl_blas.h ;;;;**************************************************************************** ;;;; Options ;;;;**************************************************************************** #+fsbv (fsbv:defcenum-aux cblas-transpose) #+fsbv (fsbv:defcenum-aux cblas-uplo) #+fsbv (fsbv:defcenum-aux cblas-diag) ;;;;**************************************************************************** ;;;; Functions ;;;;**************************************************************************** (defun matrix-product-dimensions (a b) (if (typep b 'grid:matrix) (list (first (grid:dimensions a)) (second (grid:dimensions b))) (first (grid:dimensions a)))) (defmfun matrix-product ((A grid:matrix) (x vector) &optional y (alpha 1) (beta 1) (TransA :notrans) TransB &aux (yarr (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)) :definition :generic :element-types #+fsbv :float-complex #-fsbv :float :inputs (A x) :outputs (yarr) :documentation ; FDL "If the second and third arguments are vectors, compute the matrix-vector product and sum y = alpha op(A) x + beta y, where op(A) = A, A^T, A^H 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.") (defmfun matrix-product-triangular ((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 :inputs (A x) :outputs (x) :documentation ; FDL "If the second argument is a vector, compute the matrix-vector product x = op(A) x 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 if Side is :Left and B = alpha B op(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 inverse-matrix-product ((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)) :definition :generic :element-types :float :inputs (A x yarr) :outputs (yarr) :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 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.") #+fsbv (defmfun matrix-product-hermitian ((A grid:matrix) (x vector) &optional (y (grid:make-foreign-array element-type :dimensions (matrix-product-dimensions A x) :initial-element 0)) (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)) ("gsl_blas_" :type "ger" :suffix) ((alpha :element-c-type) ((mpointer x) :pointer) ((mpointer y) :pointer) ((mpointer A) :pointer)) :definition :generic :element-types #+fsbv :float-complex #-fsbv :float :inputs (x y A) :outputs (A) :documentation ; FDL "The rank-1 update A = alpha x y^T + A of the matrix A.") #+fsbv (defmfun conjugate-rank-1-update (alpha (x vector) (y vector) (A grid:matrix)) ("gsl_blas_" :type "gerc") ((alpha :element-c-type) ((mpointer x) :pointer) ((mpointer y) :pointer) ((mpointer A) :pointer)) :definition :generic :element-types :complex :inputs (x y A) :outputs (A) :documentation ; FDL "The conjugate rank-1 update A = alpha x y^H + A of the matrix A.") (defmfun symmetric-rank-1-update ((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.") #+fsbv (defmfun hermitian-rank-1-update ((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.") (defmfun symmetric-rank-2-update ((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.") #+fsbv (defmfun hermitian-rank-2-update ((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))) (v1 (array-default 3)) (v3 (array-default 3)) (s1 (scalar-default)) (s2 (scalar-default))) (grid:copy-to (matrix-product-symmetric m1 v1 v3 s1 s2)))) #+fsbv (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))))