From 44861dcbfeee041223c4aac1ee075e92fa4daa01 Mon Sep 17 00:00:00 2001 From: Stanislaw Halik Date: Sun, 18 Sep 2016 12:42:15 +0200 Subject: update --- eigen/blas/level3_impl.h | 634 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 634 insertions(+) create mode 100644 eigen/blas/level3_impl.h (limited to 'eigen/blas/level3_impl.h') diff --git a/eigen/blas/level3_impl.h b/eigen/blas/level3_impl.h new file mode 100644 index 0000000..07dbc22 --- /dev/null +++ b/eigen/blas/level3_impl.h @@ -0,0 +1,634 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009-2010 Gael Guennebaud +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "common.h" + +int EIGEN_BLAS_FUNC(gemm)(char *opa, char *opb, int *m, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc) +{ +// std::cerr << "in gemm " << *opa << " " << *opb << " " << *m << " " << *n << " " << *k << " " << *lda << " " << *ldb << " " << *ldc << " " << *palpha << " " << *pbeta << "\n"; + typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar, internal::level3_blocking&, Eigen::internal::GemmParallelInfo*); + static functype func[12]; + + static bool init = false; + if(!init) + { + for(int k=0; k<12; ++k) + func[k] = 0; + func[NOTR | (NOTR << 2)] = (internal::general_matrix_matrix_product::run); + func[TR | (NOTR << 2)] = (internal::general_matrix_matrix_product::run); + func[ADJ | (NOTR << 2)] = (internal::general_matrix_matrix_product::run); + func[NOTR | (TR << 2)] = (internal::general_matrix_matrix_product::run); + func[TR | (TR << 2)] = (internal::general_matrix_matrix_product::run); + func[ADJ | (TR << 2)] = (internal::general_matrix_matrix_product::run); + func[NOTR | (ADJ << 2)] = (internal::general_matrix_matrix_product::run); + func[TR | (ADJ << 2)] = (internal::general_matrix_matrix_product::run); + func[ADJ | (ADJ << 2)] = (internal::general_matrix_matrix_product::run); + init = true; + } + + 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); + + int info = 0; + if(OP(*opa)==INVALID) info = 1; + else if(OP(*opb)==INVALID) info = 2; + else if(*m<0) info = 3; + else if(*n<0) info = 4; + else if(*k<0) info = 5; + else if(*lda blocking(*m,*n,*k); + + int code = OP(*opa) | (OP(*opb) << 2); + func[code](*m, *n, *k, a, *lda, b, *ldb, c, *ldc, alpha, blocking, 0); + return 0; +} + +int EIGEN_BLAS_FUNC(trsm)(char *side, char *uplo, char *opa, char *diag, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb) +{ +// std::cerr << "in trsm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << "," << *n << " " << *palpha << " " << *lda << " " << *ldb<< "\n"; + typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, internal::level3_blocking&); + static functype func[32]; + + static bool init = false; + if(!init) + { + for(int k=0; k<32; ++k) + func[k] = 0; + + func[NOTR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix::run); + func[TR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix::run); + func[ADJ | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix::run); + + func[NOTR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix::run); + func[TR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix::run); + func[ADJ | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix::run); + + func[NOTR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix::run); + func[TR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix::run); + func[ADJ | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix::run); + + func[NOTR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix::run); + func[TR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix::run); + func[ADJ | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix::run); + + + func[NOTR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix::run); + func[TR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix::run); + func[ADJ | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix::run); + + func[NOTR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix::run); + func[TR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix::run); + func[ADJ | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix::run); + + func[NOTR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix::run); + func[TR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix::run); + func[ADJ | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix::run); + + func[NOTR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix::run); + func[TR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix::run); + func[ADJ | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix::run); + + init = true; + } + + Scalar* a = reinterpret_cast(pa); + Scalar* b = reinterpret_cast(pb); + Scalar alpha = *reinterpret_cast(palpha); + + int info = 0; + if(SIDE(*side)==INVALID) info = 1; + else if(UPLO(*uplo)==INVALID) info = 2; + else if(OP(*opa)==INVALID) info = 3; + else if(DIAG(*diag)==INVALID) info = 4; + else if(*m<0) info = 5; + else if(*n<0) info = 6; + else if(*lda blocking(*m,*n,*m); + func[code](*m, *n, a, *lda, b, *ldb, blocking); + } + else + { + internal::gemm_blocking_space blocking(*m,*n,*n); + func[code](*n, *m, a, *lda, b, *ldb, blocking); + } + + if(alpha!=Scalar(1)) + matrix(b,*m,*n,*ldb) *= alpha; + + return 0; +} + + +// b = alpha*op(a)*b for side = 'L'or'l' +// b = alpha*b*op(a) for side = 'R'or'r' +int EIGEN_BLAS_FUNC(trmm)(char *side, char *uplo, char *opa, char *diag, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb) +{ +// std::cerr << "in trmm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << " " << *n << " " << *lda << " " << *ldb << " " << *palpha << "\n"; + typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&, internal::level3_blocking&); + static functype func[32]; + static bool init = false; + if(!init) + { + for(int k=0; k<32; ++k) + func[k] = 0; + + func[NOTR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + func[TR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + func[ADJ | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + + func[NOTR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + func[TR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + func[ADJ | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + + func[NOTR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + func[TR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + func[ADJ | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + + func[NOTR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + func[TR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + func[ADJ | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + + func[NOTR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + func[TR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + func[ADJ | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + + func[NOTR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + func[TR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + func[ADJ | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + + func[NOTR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + func[TR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + func[ADJ | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + + func[NOTR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + func[TR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + func[ADJ | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix::run); + + init = true; + } + + Scalar* a = reinterpret_cast(pa); + Scalar* b = reinterpret_cast(pb); + Scalar alpha = *reinterpret_cast(palpha); + + int info = 0; + if(SIDE(*side)==INVALID) info = 1; + else if(UPLO(*uplo)==INVALID) info = 2; + else if(OP(*opa)==INVALID) info = 3; + else if(DIAG(*diag)==INVALID) info = 4; + else if(*m<0) info = 5; + else if(*n<0) info = 6; + else if(*lda tmp = matrix(b,*m,*n,*ldb); + matrix(b,*m,*n,*ldb).setZero(); + + if(SIDE(*side)==LEFT) + { + internal::gemm_blocking_space blocking(*m,*n,*m); + func[code](*m, *n, *m, a, *lda, tmp.data(), tmp.outerStride(), b, *ldb, alpha, blocking); + } + else + { + internal::gemm_blocking_space blocking(*m,*n,*n); + func[code](*m, *n, *n, tmp.data(), tmp.outerStride(), a, *lda, b, *ldb, alpha, blocking); + } + return 1; +} + +// c = alpha*a*b + beta*c for side = 'L'or'l' +// c = alpha*b*a + beta*c for side = 'R'or'r +int EIGEN_BLAS_FUNC(symm)(char *side, char *uplo, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc) +{ +// std::cerr << "in symm " << *side << " " << *uplo << " " << *m << "x" << *n << " lda:" << *lda << " ldb:" << *ldb << " ldc:" << *ldc << " alpha:" << *palpha << " beta:" << *pbeta << "\n"; + 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); + + int info = 0; + if(SIDE(*side)==INVALID) info = 1; + else if(UPLO(*uplo)==INVALID) info = 2; + else if(*m<0) info = 3; + else if(*n<0) info = 4; + else if(*lda matA(size,size); + if(UPLO(*uplo)==UP) + { + matA.triangularView() = matrix(a,size,size,*lda); + matA.triangularView() = matrix(a,size,size,*lda).transpose(); + } + else if(UPLO(*uplo)==LO) + { + matA.triangularView() = matrix(a,size,size,*lda); + matA.triangularView() = matrix(a,size,size,*lda).transpose(); + } + if(SIDE(*side)==LEFT) + matrix(c, *m, *n, *ldc) += alpha * matA * matrix(b, *m, *n, *ldb); + else if(SIDE(*side)==RIGHT) + matrix(c, *m, *n, *ldc) += alpha * matrix(b, *m, *n, *ldb) * matA; + #else + if(SIDE(*side)==LEFT) + if(UPLO(*uplo)==UP) internal::product_selfadjoint_matrix::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha); + else if(UPLO(*uplo)==LO) internal::product_selfadjoint_matrix::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha); + else return 0; + else if(SIDE(*side)==RIGHT) + if(UPLO(*uplo)==UP) internal::product_selfadjoint_matrix::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha); + else if(UPLO(*uplo)==LO) internal::product_selfadjoint_matrix::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha); + else return 0; + else + return 0; + #endif + + return 0; +} + +// c = alpha*a*a' + beta*c for op = 'N'or'n' +// c = alpha*a'*a + beta*c for op = 'T'or't','C'or'c' +int EIGEN_BLAS_FUNC(syrk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc) +{ +// std::cerr << "in syrk " << *uplo << " " << *op << " " << *n << " " << *k << " " << *palpha << " " << *lda << " " << *pbeta << " " << *ldc << "\n"; + #if !ISCOMPLEX + typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&); + static functype func[8]; + + static bool init = false; + if(!init) + { + for(int k=0; k<8; ++k) + func[k] = 0; + + func[NOTR | (UP << 2)] = (internal::general_matrix_matrix_triangular_product::run); + func[TR | (UP << 2)] = (internal::general_matrix_matrix_triangular_product::run); + func[ADJ | (UP << 2)] = (internal::general_matrix_matrix_triangular_product::run); + + func[NOTR | (LO << 2)] = (internal::general_matrix_matrix_triangular_product::run); + func[TR | (LO << 2)] = (internal::general_matrix_matrix_triangular_product::run); + func[ADJ | (LO << 2)] = (internal::general_matrix_matrix_triangular_product::run); + + init = true; + } + #endif + + Scalar* a = reinterpret_cast(pa); + Scalar* c = reinterpret_cast(pc); + Scalar alpha = *reinterpret_cast(palpha); + Scalar beta = *reinterpret_cast(pbeta); + + int info = 0; + if(UPLO(*uplo)==INVALID) info = 1; + else if(OP(*op)==INVALID) info = 2; + else if(*n<0) info = 3; + else if(*k<0) info = 4; + else if(*lda().setZero(); + else matrix(c, *n, *n, *ldc).triangularView() *= beta; + else + if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView().setZero(); + else matrix(c, *n, *n, *ldc).triangularView() *= beta; + } + + #if ISCOMPLEX + // FIXME add support for symmetric complex matrix + if(UPLO(*uplo)==UP) + { + if(OP(*op)==NOTR) + matrix(c, *n, *n, *ldc).triangularView() += alpha * matrix(a,*n,*k,*lda) * matrix(a,*n,*k,*lda).transpose(); + else + matrix(c, *n, *n, *ldc).triangularView() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda); + } + else + { + if(OP(*op)==NOTR) + matrix(c, *n, *n, *ldc).triangularView() += alpha * matrix(a,*n,*k,*lda) * matrix(a,*n,*k,*lda).transpose(); + else + matrix(c, *n, *n, *ldc).triangularView() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda); + } + #else + int code = OP(*op) | (UPLO(*uplo) << 2); + func[code](*n, *k, a, *lda, a, *lda, c, *ldc, alpha); + #endif + + return 0; +} + +// c = alpha*a*b' + alpha*b*a' + beta*c for op = 'N'or'n' +// c = alpha*a'*b + alpha*b'*a + beta*c for op = 'T'or't' +int EIGEN_BLAS_FUNC(syr2k)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc) +{ + 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); + + int info = 0; + if(UPLO(*uplo)==INVALID) info = 1; + else if(OP(*op)==INVALID) info = 2; + else if(*n<0) info = 3; + else if(*k<0) info = 4; + else if(*lda().setZero(); + else matrix(c, *n, *n, *ldc).triangularView() *= beta; + else + if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView().setZero(); + else matrix(c, *n, *n, *ldc).triangularView() *= beta; + } + + if(*k==0) + return 1; + + if(OP(*op)==NOTR) + { + if(UPLO(*uplo)==UP) + { + matrix(c, *n, *n, *ldc).triangularView() + += alpha *matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).transpose() + + alpha*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).transpose(); + } + else if(UPLO(*uplo)==LO) + matrix(c, *n, *n, *ldc).triangularView() + += alpha*matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).transpose() + + alpha*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).transpose(); + } + else if(OP(*op)==TR || OP(*op)==ADJ) + { + if(UPLO(*uplo)==UP) + matrix(c, *n, *n, *ldc).triangularView() + += alpha*matrix(a, *k, *n, *lda).transpose()*matrix(b, *k, *n, *ldb) + + alpha*matrix(b, *k, *n, *ldb).transpose()*matrix(a, *k, *n, *lda); + else if(UPLO(*uplo)==LO) + matrix(c, *n, *n, *ldc).triangularView() + += alpha*matrix(a, *k, *n, *lda).transpose()*matrix(b, *k, *n, *ldb) + + alpha*matrix(b, *k, *n, *ldb).transpose()*matrix(a, *k, *n, *lda); + } + + return 0; +} + + +#if ISCOMPLEX + +// c = alpha*a*b + beta*c for side = 'L'or'l' +// c = alpha*b*a + beta*c for side = 'R'or'r +int EIGEN_BLAS_FUNC(hemm)(char *side, char *uplo, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc) +{ + 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); + +// std::cerr << "in hemm " << *side << " " << *uplo << " " << *m << " " << *n << " " << alpha << " " << *lda << " " << beta << " " << *ldc << "\n"; + + int info = 0; + if(SIDE(*side)==INVALID) info = 1; + else if(UPLO(*uplo)==INVALID) info = 2; + else if(*m<0) info = 3; + else if(*n<0) info = 4; + else if(*lda + ::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha); + else if(UPLO(*uplo)==LO) internal::product_selfadjoint_matrix + ::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha); + else return 0; + } + else if(SIDE(*side)==RIGHT) + { + if(UPLO(*uplo)==UP) matrix(c,*m,*n,*ldc) += alpha * matrix(b,*m,*n,*ldb) * matrix(a,*n,*n,*lda).selfadjointView();/*internal::product_selfadjoint_matrix + ::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);*/ + else if(UPLO(*uplo)==LO) internal::product_selfadjoint_matrix + ::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha); + else return 0; + } + else + { + return 0; + } + + return 0; +} + +// c = alpha*a*conj(a') + beta*c for op = 'N'or'n' +// c = alpha*conj(a')*a + beta*c for op = 'C'or'c' +int EIGEN_BLAS_FUNC(herk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc) +{ + typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&); + static functype func[8]; + + static bool init = false; + if(!init) + { + for(int k=0; k<8; ++k) + func[k] = 0; + + func[NOTR | (UP << 2)] = (internal::general_matrix_matrix_triangular_product::run); + func[ADJ | (UP << 2)] = (internal::general_matrix_matrix_triangular_product::run); + + func[NOTR | (LO << 2)] = (internal::general_matrix_matrix_triangular_product::run); + func[ADJ | (LO << 2)] = (internal::general_matrix_matrix_triangular_product::run); + + init = true; + } + + Scalar* a = reinterpret_cast(pa); + Scalar* c = reinterpret_cast(pc); + RealScalar alpha = *palpha; + RealScalar beta = *pbeta; + +// std::cerr << "in herk " << *uplo << " " << *op << " " << *n << " " << *k << " " << alpha << " " << *lda << " " << beta << " " << *ldc << "\n"; + + int info = 0; + if(UPLO(*uplo)==INVALID) info = 1; + else if((OP(*op)==INVALID) || (OP(*op)==TR)) info = 2; + else if(*n<0) info = 3; + else if(*k<0) info = 4; + else if(*lda().setZero(); + else matrix(c, *n, *n, *ldc).triangularView() *= beta; + else + if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView().setZero(); + else matrix(c, *n, *n, *ldc).triangularView() *= beta; + + if(beta!=Scalar(0)) + { + matrix(c, *n, *n, *ldc).diagonal().real() *= beta; + matrix(c, *n, *n, *ldc).diagonal().imag().setZero(); + } + } + + if(*k>0 && alpha!=RealScalar(0)) + { + func[code](*n, *k, a, *lda, a, *lda, c, *ldc, alpha); + matrix(c, *n, *n, *ldc).diagonal().imag().setZero(); + } + return 0; +} + +// c = alpha*a*conj(b') + conj(alpha)*b*conj(a') + beta*c, for op = 'N'or'n' +// c = alpha*conj(a')*b + conj(alpha)*conj(b')*a + beta*c, for op = 'C'or'c' +int EIGEN_BLAS_FUNC(her2k)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc) +{ + Scalar* a = reinterpret_cast(pa); + Scalar* b = reinterpret_cast(pb); + Scalar* c = reinterpret_cast(pc); + Scalar alpha = *reinterpret_cast(palpha); + RealScalar beta = *pbeta; + + int info = 0; + if(UPLO(*uplo)==INVALID) info = 1; + else if((OP(*op)==INVALID) || (OP(*op)==TR)) info = 2; + else if(*n<0) info = 3; + else if(*k<0) info = 4; + else if(*lda().setZero(); + else matrix(c, *n, *n, *ldc).triangularView() *= beta; + else + if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView().setZero(); + else matrix(c, *n, *n, *ldc).triangularView() *= beta; + + if(beta!=Scalar(0)) + { + matrix(c, *n, *n, *ldc).diagonal().real() *= beta; + matrix(c, *n, *n, *ldc).diagonal().imag().setZero(); + } + } + else if(*k>0 && alpha!=Scalar(0)) + matrix(c, *n, *n, *ldc).diagonal().imag().setZero(); + + if(*k==0) + return 1; + + if(OP(*op)==NOTR) + { + if(UPLO(*uplo)==UP) + { + matrix(c, *n, *n, *ldc).triangularView() + += alpha *matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint() + + numext::conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint(); + } + else if(UPLO(*uplo)==LO) + matrix(c, *n, *n, *ldc).triangularView() + += alpha*matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint() + + numext::conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint(); + } + else if(OP(*op)==ADJ) + { + if(UPLO(*uplo)==UP) + matrix(c, *n, *n, *ldc).triangularView() + += alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb) + + numext::conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda); + else if(UPLO(*uplo)==LO) + matrix(c, *n, *n, *ldc).triangularView() + += alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb) + + numext::conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda); + } + + return 1; +} + +#endif // ISCOMPLEX -- cgit v1.2.3