Source code for pygsl.blas

# Author : Fabian Jakobs
r"""BLAS - Basic Linear Algebra Subprograms

GSL provides dense vector and matrix objects, based on the relevant built-in
types. The library provides an interface to the BLAS operations which apply to
these objects.

PyGSL only provides the functions working on "native" Python datatypes, i.e.
double and complex_double.

Functions with the postfix "_cr" are functions that support call by reference.
This is faster than the original version but may confuse real Python
programmers. Use this, if speed matters!!

Functions that are naturally done using functions of the underlying numerical
package are left out here.
"""

from . import gslwrap
from . import _gslwrap
import pygsl._numobj as Numeric
import copy
from pygsl import array_typed_copy, get_typecode, Float, Complex
#
# constants
#
CblasRowMajor = gslwrap.CblasRowMajor
CblasColMajor = gslwrap.CblasColMajor
CblasNoTrans = gslwrap.CblasNoTrans
CblasTrans = gslwrap.CblasTrans
CblasConjTrans = gslwrap.CblasConjTrans
CblasUpper = gslwrap.CblasUpper
CblasLower = gslwrap.CblasLower
CblasNonUnit = gslwrap.CblasNonUnit
CblasUnit = gslwrap.CblasUnit
CblasLeft = gslwrap.CblasLeft
CblasRight = gslwrap.CblasRight
gsl_blas_sdsdot = gslwrap.gsl_blas_sdsdot

#
# Level 1
#

[docs] def ddot(x, y): r"""This function computes the scalar product \M{x^T y} for the vectors x and y, returning the result. """ return _gslwrap.gsl_blas_ddot(x, y)#[1]
[docs] def zdotu(x, y): r"""This function computes the complex scalar product \M{x^T y} for the vectors x and y, returning the result. """ return _gslwrap.gsl_blas_zdotu(x, y, 1j)#[1]
[docs] def zdotc(x, y): r"""This function computes the complex conjugate scalar product \M{x^H y} for the vectors x and y, returning the result. """ return _gslwrap.gsl_blas_zdotc(x, y, 1j)#[1]
[docs] def dnrm2(x): r"""This function computes the Euclidean norm \M{||x||_2 = \sqrt {\sum x_i^2}} of the vector x. """ return _gslwrap.gsl_blas_dnrm2(x)
[docs] def dznrm2(x): r"""This function computes the Euclidean norm of the complex vector x, \M{||x||_2 = \sqrt {\sum (\Re(x_i)^2 + \Im(x_i)^2)}}. """ return _gslwrap.gsl_blas_dznrm2(x)
[docs] def dasum(x): r"""This function computes the absolute sum \M{\sum |x_i|} of the elements of the vector x. """ return _gslwrap.gsl_blas_dasum(x)
[docs] def dzasum(x): r"""This function computes the absolute sum \M{\sum |\Re(x_i)| + |\Im(x_i)|} of the elements of the vector x. """ return _gslwrap.gsl_blas_dzasum(x)
[docs] def idamax(x): """This function returns the index of the largest element of the vector x. The largest element is determined by its absolute magnitude. If the largest value occurs several times then the index of the first occurrence is returned. """ return _gslwrap.gsl_blas_idamax(x)
[docs] def izamax(x): r"""This function returns the index of the largest element of the vector x. The largest element is determined by the sum of the magnitudes of the real and imaginary parts \M{|\Re(x_i)| + |\Im(x_i)|}. If the largest value occurs several times then the index of the first occurrence is returned. """ return _gslwrap.gsl_blas_izamax(x)
[docs] def daxpy(alpha, x, y): r"""This function computes the sum \M{y = S{alpha} x + y} for the vectors x and y. """ yn = array_typed_copy(y, Float) _gslwrap.gsl_blas_daxpy(alpha, x, yn) return yn
[docs] def daxpy_cr(alpha, x, y_CR): r"""This function computes the sum \M{y = S{alpha} x + y} for the vectors x and y. """ _gslwrap.gsl_blas_daxpy(alpha, x, y_CR)
[docs] def zaxpy(alpha, x, y): r"""This function computes the sum \M{y = S{alpha} x + y} for the vectors x and y. """ yn = array_typed_copy(y, Complex) _gslwrap.gsl_blas_zaxpy(alpha, x, yn) return yn
[docs] def zaxpy_cr(alpha, x, y_CR): r"""This function computes the sum \M{y = S{alpha} x + y} for the vectors x and y. """ _gslwrap.gsl_blas_zaxpy(alpha, x, y_CR)
[docs] def drot(x, y, c, s): r"""This function applies a Givens rotation \M{(x', y') = (c x + s y, -s x + c y)} to the vectors x, y. """ xn = array_typed_copy(x, Float) yn = array_typed_copy(y, Float) _gslwrap.gsl_blas_drot(xn, yn, c, s) return (xn, yn)
[docs] def drot_cr(x_CR, y_CR, c, s): r"""This function applies a Givens rotation \M{(x', y') = (c x + s y, -s x + c y)} to the vectors x, y. """ _gslwrap.gsl_blas_drot(x_CR, y_CR, c, s)
# # Level 2 #
[docs] def dgemv(alpha, a, x, beta, y, TransA=CblasNoTrans): r"""This function computes the matrix-vector product and sum \M{y = S{alpha} op(A) x + S{beta} y}, where op(A) = \M{A, A^T, A^H} for TransA = CblasNoTrans, CblasTrans, CblasConjTrans. """ yn = array_typed_copy(y, Float) _gslwrap.gsl_blas_dgemv(TransA, alpha, a, x, beta, yn) return yn
[docs] def zgemv(alpha, a, x, beta, y, TransA=CblasNoTrans): r"""This function computes the matrix-vector product and sum \M{y = S{alpha} op(A) x + S{beta} y}, where \M{op(A) = A, A^T, A^H} for TransA = CblasNoTrans, CblasTrans, CblasConjTrans. """ yn = array_typed_copy(y, Complex) _gslwrap.gsl_blas_zgemv(TransA, alpha, a, x, beta, yn) return yn
[docs] def dtrmv(A, x, Uplo=CblasLower, TransA=CblasNoTrans, Diag=CblasNonUnit): """This function computes the matrix-vector product returns x' This function computes the matrix-vector product and x' = op(A) x for the triangular matrix A, where op(A) = A, A^T, A^H for TransA = CblasNoTrans, CblasTrans, CblasConjTrans. When Uplo is CblasUpper then the upper triangle of A is used, and when Uplo is CblasLower then the lower triangle of A is used. If Diag is CblasNonUnit then the diagonal of the matrix is used, but if Diag is CblasUnit then the diagonal elements of the matrix A are taken as unity and are not referenced. """ xn = array_typed_copy(x, Float) _gslwrap.gsl_blas_dtrmv(Uplo, TransA, Diag, A, xn) return xn
[docs] def ztrmv(A, x, Uplo=CblasLower, TransA=CblasNoTrans, Diag=CblasNonUnit): """ returns x' This function computes the matrix-vector product and x' = op(A) x for the triangular matrix A, where op(A) = A, A^T, A^H for TransA = CblasNoTrans, CblasTrans, CblasConjTrans. When Uplo is CblasUpper then the upper triangle of A is used, and when Uplo is CblasLower then the lower triangle of A is used. If Diag is CblasNonUnit then the diagonal of the matrix is used, but if Diag is CblasUnit then the diagonal elements of the matrix A are taken as unity and are not referenced. """ xn = array_typed_copy(x) _gslwrap.gsl_blas_ztrmv(Uplo, TransA, Diag, A, xn) return xn
[docs] def dtrsv(A, x, Uplo=CblasLower, TransA=CblasNoTrans, Diag=CblasNonUnit): """ This function computes inv(op(A)) x for x, where op(A) = A, A^T, A^H for TransA = CblasNoTrans, CblasTrans, CblasConjTrans. When Uplo is CblasUpper then the upper triangle of A is used, and when Uplo is CblasLower then the lower triangle of A is used. If Diag is CblasNonUnit then the diagonal of the matrix is used, but if Diag is CblasUnit then the diagonal elements of the matrix A are taken as unity and are not referenced. Returns: x: """ xn = array_typed_copy(x) _gslwrap.gsl_blas_dtrsv(Uplo, TransA, Diag, A, xn) return xn
[docs] def ztrsv(A, x, Uplo=CblasLower, TransA=CblasNoTrans, Diag=CblasNonUnit): """ returns x' This function computes inv(op(A)) x for x, where op(A) = A, A^T, A^H for TransA = CblasNoTrans, CblasTrans, CblasConjTrans. When Uplo is CblasUpper then the upper triangle of A is used, and when Uplo is CblasLower then the lower triangle of A is used. If Diag is CblasNonUnit then the diagonal of the matrix is used, but if Diag is CblasUnit then the diagonal elements of the matrix A are taken as unity and are not referenced. """ xn = array_typed_copy(x) _gslwrap.gsl_blas_ztrsv(Uplo, TransA, Diag, A, xn) return xn
[docs] def dsymv(alpha, A, X, beta, Y, Uplo=CblasLower): r""" returns y' This function computes the matrix-vector product and sum \M{y' = S{alpha} A x + S{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 CblasUpper then the upper triangle and diagonal of A are used, and when Uplo is CblasLower then the lower triangle and diagonal of A are used. """ yn = array_typed_copy(Y, Float) _gslwrap.gsl_blas_dsymv(Uplo, alpha, A, X, beta, yn) return yn
[docs] def zhemv(alpha, A, X, beta, Y, Uplo=CblasLower): r""" returns y' This function computes the matrix-vector product and sum \M{y' = S{alpha} A x + S{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 CblasUpper then the upper triangle and diagonal of A are used, and when Uplo is CblasLower 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. """ yn = array_typed_copy(Y, get_typecode(A)) _gslwrap.gsl_blas_zhemv(Uplo, alpha, A, X, beta, yn) return yn
[docs] def dger(alpha, X, Y, A): r""" returns A' This function computes the rank-1 update \M{A' = S{alpha} x y^T + A} of the matrix A. """ an = array_typed_copy(A) _gslwrap.gsl_blas_dger(alpha, X, Y, an) return an
[docs] def zgeru(alpha, X, Y, A): r""" returns A' This function computes the rank-1 update \M{A' = S{alpha} x y^T + A} of the matrix A. """ an = array_typed_copy(A) _gslwrap.gsl_blas_zgeru(alpha, X, Y, an) return an
[docs] def zgerc(alpha, X, Y, A): r""" returns A' This function computes the conjugate rank-1 update \M{A = S{alpha} x y^H + A} of the matrix A. """ an = array_typed_copy(A) _gslwrap.gsl_blas_zgerc(alpha, X, Y, an) return an
[docs] def dsyr(alpha, X, A, Uplo=CblasLower): r""" returns A' This function computes the symmetric rank-1 update \M{A' = S{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 CblasUpper then the upper triangle and diagonal of A are used, and when Uplo is CblasLower then the lower triangle and diagonal of A are used. """ an = array_typed_copy(A) _gslwrap.gsl_blas_dsyr(Uplo, alpha, X, an) return an
[docs] def zher(alpha, X, A, Uplo=CblasLower): r""" returns A' This function computes the hermitian rank-1 update \M{A' = S{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 CblasUpper then the upper triangle and diagonal of A are used, and when Uplo is CblasLower then the lower triangle and diagonal of A are used. The imaginary elements of the diagonal are automatically set to zero. """ an = array_typed_copy(A) _gslwrap.gsl_blas_zher(Uplo, alpha, X, an) return an
[docs] def dsyr2(alpha, X, Y, A, Uplo=CblasLower): r""" returns A' This function computes the symmetric rank-2 update \M{A' = S{alpha} x y^T + S{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 CblasUpper then the upper triangle and diagonal of A are used, and when Uplo is CblasLower then the lower triangle and diagonal of A are used. """ an = array_typed_copy(A) _gslwrap.gsl_blas_dsyr2(Uplo, alpha, X, Y, an) return an
[docs] def zher2(alpha, X, Y, A, Uplo=CblasLower): r""" returns A' This function computes the hermitian rank-2 update \M{A' = S{alpha} x y^H + S{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 CblasUpper then the upper triangle and diagonal of A are used, and when Uplo is CblasLower then the lower triangle and diagonal of A are used. The imaginary elements of the diagonal are automatically set to zero. """ an = array_typed_copy(A) _gslwrap.gsl_blas_zher2(Uplo, alpha, X, Y, an) return an
# # Level 3 #
[docs] def dgemm(alpha, A, B, beta, C, TransA=CblasNoTrans, TransB=CblasNoTrans): r""" returns C' This function computes the matrix-matrix product and sum \M{C' = S{alpha} op(A) op(B) + S{beta} C} where op(A) = A, A^T, A^H for TransA = CblasNoTrans, CblasTrans, CblasConjTrans and similarly for the parameter TransB. """ cn = array_typed_copy(C) _gslwrap.gsl_blas_dgemm(TransA, TransB, alpha, A, B, beta, cn) return cn
[docs] def zgemm(alpha, A, B, beta, C, TransA=CblasNoTrans, TransB=CblasNoTrans): r""" returns C' This function computes the matrix-matrix product and sum \M{C' = S{alpha} op(A) op(B) + S{beta} C} where op(A) = A, A^T, A^H for TransA = CblasNoTrans, CblasTrans, CblasConjTrans and similarly for the parameter TransB. """ cn = array_typed_copy(C) _gslwrap.gsl_blas_zgemm(TransA, TransB, alpha, A, B, beta, cn) return cn
[docs] def dsymm(alpha, A, B, beta, C, Side=CblasLeft, Uplo=CblasLower): r""" returns C' This function computes the matrix-matrix product and sum \M{C' = S{alpha} A B + S{beta} C} for Side is CblasLeft and \M{C' = S{alpha} B A + S{beta} C} for Side is CblasRight, where the matrix A is symmetric. When Uplo is CblasUpper then the upper triangle and diagonal of A are used, and when Uplo is CblasLower then the lower triangle and diagonal of A are used. """ cn = array_typed_copy(C) _gslwrap.gsl_blas_dsymm(Side, Uplo, alpha, A, B, beta, cn) return cn
[docs] def zsymm(alpha, A, B, beta, C, Side=CblasLeft, Uplo=CblasLower): r""" returns C' This function computes the matrix-matrix product and sum \M{C' = S{alpha} A B + S{beta} C} for Side is CblasLeft and \M{C' = S{alpha} B A + S{beta} C} for Side is CblasRight, where the matrix A is symmetric. When Uplo is CblasUpper then the upper triangle and diagonal of A are used, and when Uplo is CblasLower then the lower triangle and diagonal of A are used. """ cn = array_typed_copy(C) _gslwrap.gsl_blas_zsymm(Side, Uplo, alpha, A, B, beta, cn) return cn
[docs] def zhemm(alpha, A, B, beta, C, Side=CblasLeft, Uplo=CblasLower): r""" returns C' This function computes the matrix-matrix product and sum \M{C' = S{alpha} A B + S{beta} C} for Side is CblasLeft and \M{C' = S{alpha} B A + S{beta} C} for Side is CblasRight, where the matrix A is hermitian. When Uplo is CblasUpper then the upper triangle and diagonal of A are used, and when Uplo is CblasLower then the lower triangle and diagonal of A are used. The imaginary elements of the diagonal are automatically set to zero. """ cn = array_typed_copy(C) _gslwrap.gsl_blas_zhemm(Side, Uplo, alpha, A, B, beta, cn) return cn
[docs] def dtrmm(alpha, A, B, Side=CblasLeft, Uplo=CblasLower, TransA=CblasNoTrans, Diag=CblasNonUnit): r""" returns B' This function computes the matrix-matrix product \M{B' = S{alpha} op(A) B} for Side is CblasLeft and \M{B' = S{alpha} B op(A)} for Side is CblasRight. The matrix A is triangular and op(A) = A, A^T, A^H for TransA = CblasNoTrans, CblasTrans, CblasConjTrans When Uplo is CblasUpper then the upper triangle of A is used, and when Uplo is CblasLower then the lower triangle of A is used. If Diag is CblasNonUnit then the diagonal of A is used, but if Diag is CblasUnit then the diagonal elements of the matrix A are taken as unity and are not referenced. """ bn = array_typed_copy(B) _gslwrap.gsl_blas_dtrmm(Side, Uplo, TransA, Diag, alpha, A, bn) return bn
[docs] def ztrmm(alpha, A, B, Side=CblasLeft, Uplo=CblasLower, TransA=CblasNoTrans, Diag=CblasNonUnit): r""" returns B' This function computes the matrix-matrix product \M{B' = S{alpha} op(A) B} for Side is CblasLeft and \M{B' = S{alpha} B op(A)} for Side is CblasRight. The matrix A is triangular and op(A) = A, A^T, A^H for TransA = CblasNoTrans, CblasTrans, CblasConjTrans When Uplo is CblasUpper then the upper triangle of A is used, and when Uplo is CblasLower then the lower triangle of A is used. If Diag is CblasNonUnit then the diagonal of A is used, but if Diag is CblasUnit then the diagonal elements of the matrix A are taken as unity and are not referenced. """ bn = array_typed_copy(B) _gslwrap.gsl_blas_ztrmm(Side, Uplo, TransA, Diag, alpha, A, bn) return bn
[docs] def dtrsm(alpha, A, B, Side=CblasLeft, Uplo=CblasLower, TransA=CblasNoTrans, Diag=CblasNonUnit): r""" returns B' This function computes the matrix-matrix product \M{B' = S{alpha} op(inv(A)) B} for Side is CblasLeft and \M{B' = S{alpha} B op(inv(A))} for Side is CblasRight. The matrix A is triangular and op(A) = A, A^T, A^H for TransA = CblasNoTrans, CblasTrans, CblasConjTrans When Uplo is CblasUpper then the upper triangle of A is used, and when Uplo is CblasLower then the lower triangle of A is used. If Diag is CblasNonUnit then the diagonal of A is used, but if Diag is CblasUnit then the diagonal elements of the matrix A are taken as unity and are not referenced. """ bn = array_typed_copy(B) _gslwrap.gsl_blas_dtrsm(Side, Uplo, TransA, Diag, alpha, A, bn) return bn
[docs] def ztrsm(alpha, A, B, Side=CblasLeft, Uplo=CblasLower, TransA=CblasNoTrans, Diag=CblasNonUnit): r""" Returns: :math:`B'` This function computes the matrix-matrix product \M{B' = S{alpha} op(inv(A)) B} for Side is CblasLeft and \M{B' = S{alpha} B op(inv(A))} for Side is CblasRight. The matrix A is triangular and op(A) = A, A^T, A^H for TransA = CblasNoTrans, CblasTrans, CblasConjTrans When Uplo is CblasUpper then the upper triangle of A is used, and when Uplo is CblasLower then the lower triangle of A is used. If Diag is CblasNonUnit then the diagonal of A is used, but if Diag is CblasUnit then the diagonal elements of the matrix A are taken as unity and are not referenced. """ bn = array_typed_copy(B) _gslwrap.gsl_blas_ztrsm(Side, Uplo, TransA, Diag, alpha, A, bn) return bn
[docs] def dsyrk(alpha, A, beta, C, Uplo=CblasLower, Trans=CblasNoTrans): r""" returns C' This function computes a rank-k update of the symmetric matrix C, \M{C' = S{alpha} A A^T + S{beta} C} when Trans is CblasNoTrans and \M{C' = S{alpha} A^T A + S{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. """ cn = array_typed_copy(C) _gslwrap.gsl_blas_dsyrk(Uplo, Trans, alpha, A, beta, cn) return cn
[docs] def zsyrk(alpha, A, beta, C, Uplo=CblasLower, Trans=CblasNoTrans): r""" returns C' This function computes a rank-k update of the symmetric matrix C, \M{C' = S{alpha} A A^T + S{beta} C} when Trans is CblasNoTrans and \M{C' = S{alpha} A^T A + S{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. """ cn = array_typed_copy(C) _gslwrap.gsl_blas_zsyrk(Uplo, Trans, alpha, A, beta, cn) return cn
[docs] def triang2symm(A, Uplo=CblasLower, Diag=CblasNonUnit): """ returns A' This function creates an symmetric matrix from a triangular matrix. When Uplo is CblasUpper then the upper triangle and diagonal of A are used, and when Uplo is CblasLower then the lower triangle and diagonal of A are used. If Diag is CblasNonUnit then the diagonal of A is used, but if Diag is CblasUnit then the diagonal elements of the matrix A are taken as unity and are not referenced. """ an = array_typed_copy(A) if Uplo == CblasLower: for i in range(A.shape[0]): an[:i+1,i] = A[i,:i+1] else: for i in range(A.shape[0]): an[i:,i] = A[i,i:] if Diag == CblasUnit: for i in range(an.shape[0]): an[i,i] = 1 return an
[docs] def triang2herm(A, Uplo=CblasLower, Diag=CblasNonUnit): """ returns A' This function creates an hermitian matrix from a triangular matrix. When Uplo is CblasUpper then the upper triangle and diagonal of A are used, and when Uplo is CblasLower then the lower triangle and diagonal of A are used. If Diag is CblasNonUnit then the diagonal of A is used, but if Diag is CblasUnit then the diagonal elements of the matrix A are taken as unity and are not referenced. """ an = array_typed_copy(A) if Uplo == CblasLower: for i in range(A.shape[0]): an[:i+1,i] = Numeric.conjugate(A[i,:i+1]) else: for i in range(A.shape[0]): an[i:,i] = Numeric.conjugate(A[i,i:]) if Diag == CblasUnit: for i in range(an.shape[0]): an[i,i] = 1 return an
[docs] def zherk(alpha, A, beta, C, Uplo=CblasLower, Trans=CblasNoTrans): r""" returns C' This function computes a rank-k update of the hermitian matrix C, \M{C' = S{alpha} A A^H + S{beta} C} when Trans is CblasNoTrans and \M{C' = S{alpha} A^H A + S{beta} C} when Trans is CblasTrans. Since the matrix C is hermitian 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. The imaginary elements of the diagonal are automatically set to zero. """ cn = array_typed_copy(C) _gslwrap.gsl_blas_zherk(Uplo, Trans, alpha, A, beta, cn) return cn
[docs] def dsyr2k(alpha, A, B, beta, C, Uplo=CblasLower, Trans=CblasNoTrans): r""" returns C' This function computes a rank-2k update of the symmetric matrix C, \M{C' = S{alpha} A B^T + S{alpha} B A^T + S{beta} C} when Trans is CblasNoTrans and \M{C' = S{alpha} A^T B + S{alpha} B^T A + S{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. """ cn = array_typed_copy(C) _gslwrap.gsl_blas_dsyr2k(Uplo, Trans, alpha, A, B, beta, cn) return cn
[docs] def zsyr2k(alpha, A, B, beta, C, Uplo=CblasLower, Trans=CblasNoTrans): r""" returns C' This function computes a rank-2k update of the symmetric matrix C, \M{C' = S{alpha} A B^T + S{alpha} B A^T + S{beta} C} when Trans is CblasNoTrans and \M{C' = S{alpha} A^T B + S{alpha} B^T A + S{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. """ cn = array_typed_copy(C) _gslwrap.gsl_blas_zsyr2k(Uplo, Trans, alpha, A, B, beta, cn) return cn
[docs] def zher2k(alpha, A, B, beta, C, Uplo=CblasLower, Trans=CblasNoTrans): r""" returns C' This function computes a rank-2k update of the hermitian matrix C, \M{C' = S{alpha} A B^H + S{alpha}^* B A^H + S{beta} C} when Trans is CblasNoTrans and \M{C' = S{alpha} A^H B + S{alpha}^* B^H A + S{beta} C} when Trans is CblasTrans. Since the matrix C is hermitian 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. The imaginary elements of the diagonal are automatically set to zero. Calls function: """ cn = array_typed_copy(C) _gslwrap.gsl_blas_zher2k(Uplo, Trans, alpha, A, B, beta, cn) return cn