From 04dc63776a63e5d0ec0237706cb440152d57769e Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 25 Sep 2009 13:08:39 +0200 Subject: add a wip blas library built on top of Eigen. TODO: - write extentive unit tests (maybe this already exist in other projects) - the level2 functions still have to be implemented --- blas/level2_impl.h | 214 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 214 insertions(+) create mode 100644 blas/level2_impl.h (limited to 'blas/level2_impl.h') diff --git a/blas/level2_impl.h b/blas/level2_impl.h new file mode 100644 index 000000000..5691e8a7f --- /dev/null +++ b/blas/level2_impl.h @@ -0,0 +1,214 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, 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 2 of +// the License, or (at your option) any later version. +// +// Eigen 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 Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#include "common.h" + +int EIGEN_BLAS_FUNC(gemv)(char *opa, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *incb, RealScalar *pbeta, RealScalar *pc, int *incc) +{ + Scalar* a = reinterpret_cast(pa); + Scalar* b = reinterpret_cast(pb); + Scalar* c = reinterpret_cast(pc); + Scalar alpha = *reinterpret_cast(palpha); + Scalar beta = *reinterpret_cast(pbeta); + + if(beta!=Scalar(1)) + vector(c, *m, *incc) *= beta; + + if(OP(*opa)==NOTR) + if(*incc==1) + vector(c,*m) += alpha * matrix(a,*m,*n,*lda) * vector(b,*n,*incb); + else + vector(c,*m,*incc) += alpha * matrix(a,*m,*n,*lda) * vector(b,*n,*incb); + else if(OP(*opa)==TR) + if(*incb==1) + vector(c,*m,*incc) += alpha * matrix(a,*n,*m,*lda).transpose() * vector(b,*n); + else + vector(c,*m,*incc) += alpha * matrix(a,*n,*m,*lda).transpose() * vector(b,*n,*incb); + else if(OP(*opa)==TR) + if(*incb==1) + vector(c,*m,*incc) += alpha * matrix(a,*n,*m,*lda).adjoint() * vector(b,*n); + else + vector(c,*m,*incc) += alpha * matrix(a,*n,*m,*lda).adjoint() * vector(b,*n,*incb); + else + return 0; + + return 1; +} + +/* +int EIGEN_BLAS_FUNC(trsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb) +{ + typedef void (*functype)(int, const Scalar *, int, Scalar *, int); + functype func[16]; + + static bool init = false; + if(!init) + { + for(int k=0; k<16; ++k) + func[k] = 0; + +// func[NOTR | (UP << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector::run); +// func[TR | (UP << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector::run); +// func[ADJ | (UP << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector::run); +// +// func[NOTR | (LO << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector::run); +// func[TR | (LO << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector::run); +// func[ADJ | (LO << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector::run); +// +// func[NOTR | (UP << 3) | (UNIT << 3)] = (ei_triangular_solve_vector::run); +// func[TR | (UP << 2) | (UNIT << 3)] = (ei_triangular_solve_vector::run); +// func[ADJ | (UP << 2) | (UNIT << 3)] = (ei_triangular_solve_vector::run); +// +// func[NOTR | (LO << 2) | (UNIT << 3)] = (ei_triangular_solve_vector::run); +// func[TR | (LO << 2) | (UNIT << 3)] = (ei_triangular_solve_vector::run); +// func[ADJ | (LO << 2) | (UNIT << 3)] = (ei_triangular_solve_vector::run); + + init = true; + } + + Scalar* a = reinterpret_cast(pa); + Scalar* b = reinterpret_cast(pb); + + int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); + if(code>=16 || func[code]==0) + return 0; + + func[code](*n, a, *lda, b, *incb); + return 1; +} +*/ + +/* +int EIGEN_BLAS_FUNC(trmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb) +{ + // TODO + + typedef void (*functype)(int, const Scalar *, int, const Scalar *, int, Scalar *, int); + functype func[16]; + + static bool init = false; + if(!init) + { + for(int k=0; k<16; ++k) + func[k] = 0; + +// func[NOTR | (UP << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector::run); +// func[TR | (UP << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector::run); +// func[ADJ | (UP << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector::run); +// +// func[NOTR | (LO << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector::run); +// func[TR | (LO << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector::run); +// func[ADJ | (LO << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector::run); +// +// func[NOTR | (UP << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector::run); +// func[TR | (UP << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector::run); +// func[ADJ | (UP << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector::run); +// +// func[NOTR | (LO << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector::run); +// func[TR | (LO << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector::run); +// func[ADJ | (LO << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector::run); + + init = true; + } + + Scalar* a = reinterpret_cast(pa); + Scalar* b = reinterpret_cast(pb); + + int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); + if(code>=16 || func[code]==0) + return 0; + + func[code](*n, a, *lda, b, *incb, b, *incb); + return 1; +} +*/ + +/* +int EIGEN_BLAS_FUNC(syr)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *inca, RealScalar *pc, int *ldc) +{ + // TODO + typedef void (*functype)(int, const Scalar *, int, Scalar *, int, Scalar); + functype func[2]; + + static bool init = false; + if(!init) + { + for(int k=0; k<2; ++k) + func[k] = 0; + +// func[UP] = (ei_selfadjoint_product::run); +// func[LO] = (ei_selfadjoint_product::run); + + init = true; + } + + Scalar* a = reinterpret_cast(pa); + Scalar* c = reinterpret_cast(pc); + Scalar alpha = *reinterpret_cast(palpha); + + int code = UPLO(*uplo); + if(code>=2 || func[code]==0) + return 0; + + func[code](*n, a, *inca, c, *ldc, alpha); + return 1; +} +*/ + +/* +int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *inca, RealScalar *pb, int *incb, RealScalar *pc, int *ldc) +{ + // TODO + typedef void (*functype)(int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar); + functype func[2]; + + static bool init = false; + if(!init) + { + for(int k=0; k<2; ++k) + func[k] = 0; + +// func[UP] = (ei_selfadjoint_product::run); +// func[LO] = (ei_selfadjoint_product::run); + + init = true; + } + + Scalar* a = reinterpret_cast(pa); + Scalar* b = reinterpret_cast(pb); + Scalar* c = reinterpret_cast(pc); + Scalar alpha = *reinterpret_cast(palpha); + + int code = UPLO(*uplo); + if(code>=2 || func[code]==0) + return 0; + + func[code](*n, a, *inca, b, *incb, c, *ldc, alpha); + return 1; +} +*/ + +#if ISCOMPLEX + +#endif // ISCOMPLEX -- cgit v1.2.3