Quelle level3_impl.h Sprache: C

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// 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 <iostream>
#include "common.h"

int EIGEN_BLAS_FUNC(gemm)(const char *opa, const char *opb, const int *m, const int *n, const int *k, const RealScalar *palpha,
                          const RealScalar *pa, const int *lda, const RealScalar *pb, const int *ldb, const RealScalar *pbeta, RealScalar *pc, const 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, DenseIndex, Scalar, internal::level3_blocking<Scalar,Scalar>&, Eigen::internal::GemmParallelInfo<DenseIndex>*);
  static const functype func[12] = {
    // array index: NOTR  | (NOTR << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,ColMajor,false,ColMajor,1>::run),
    // array index: TR    | (NOTR << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,false,ColMajor,1>::run),
    // array index: ADJ   | (NOTR << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,1>::run),
    0,
    // array index: NOTR  | (TR   << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,false,ColMajor,1>::run),
    // array index: TR    | (TR   << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,RowMajor,false,ColMajor,1>::run),
    // array index: ADJ   | (TR   << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,RowMajor,false,ColMajor,1>::run),
    0,
    // array index: NOTR  | (ADJ  << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,1>::run),
    // array index: TR    | (ADJ  << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,RowMajor,Conj, ColMajor,1>::run),
    // array index: ADJ   | (ADJ  << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,RowMajor,Conj, ColMajor,1>::run),
    0
  };

  const Scalar* a = reinterpret_cast<const Scalar*>(pa);
  const Scalar* b = reinterpret_cast<const Scalar*>(pb);
  Scalar* c = reinterpret_cast<Scalar*>(pc);
  Scalar alpha  = *reinterpret_cast<const Scalar*>(palpha);
  Scalar beta   = *reinterpret_cast<const Scalar*>(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<std::max(1,(OP(*opa)==NOTR)?*m:*k))                    info = 8;
  else if(*ldb<std::max(1,(OP(*opb)==NOTR)?*k:*n))                    info = 10;
  else if(*ldc<std::max(1,*m))                                        info = 13;
  if(info)
    return xerbla_(SCALAR_SUFFIX_UP"GEMM ",&info,6);

  if (*m == 0 || *n == 0)
    return 0;

  if(beta!=Scalar(1))
  {
    if(beta==Scalar(0)) matrix(c, *m, *n, *ldc).setZero();
    else                matrix(c, *m, *n, *ldc) *= beta;
  }

  if(*k == 0)
    return 0;

  internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic> blocking(*m,*n,*k,1,true);

  int code = OP(*opa) | (OP(*opb) << 2);
  func[code](*m, *n, *k, a, *lda, b, *ldb, c, 1, *ldc, alpha, blocking, 0);
  return 0;
}

int EIGEN_BLAS_FUNC(trsm)(const char *side, const char *uplo, const char *opa, const char *diag, const int *m, const int *n,
                          const RealScalar *palpha,  const RealScalar *pa, const int *lda, RealScalar *pb, const 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, DenseIndex, internal::level3_blocking<Scalar,Scalar>&);
  static const functype func[32] = {
    // array index: NOTR  | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0,          false,ColMajor,ColMajor,1>::run),
    // array index: TR    | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0,          false,RowMajor,ColMajor,1>::run),
    // array index: ADJ   | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0,          Conj, RowMajor,ColMajor,1>::run),\
    0,
    // array index: NOTR  | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0,          false,ColMajor,ColMajor,1>::run),
    // array index: TR    | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0,          false,RowMajor,ColMajor,1>::run),
    // array index: ADJ   | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0,          Conj, RowMajor,ColMajor,1>::run),
    0,
    // array index: NOTR  | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0,          false,ColMajor,ColMajor,1>::run),
    // array index: TR    | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0,          false,RowMajor,ColMajor,1>::run),
    // array index: ADJ   | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0,          Conj, RowMajor,ColMajor,1>::run),
    0,
    // array index: NOTR  | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0,          false,ColMajor,ColMajor,1>::run),
    // array index: TR    | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0,          false,RowMajor,ColMajor,1>::run),
    // array index: ADJ   | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0,          Conj, RowMajor,ColMajor,1>::run),
    0,
    // array index: NOTR  | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,false,ColMajor,ColMajor,1>::run),
    // array index: TR    | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,false,RowMajor,ColMajor,1>::run),
    // array index: ADJ   | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,Conj, RowMajor,ColMajor,1>::run),
    0,
    // array index: NOTR  | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,false,ColMajor,ColMajor,1>::run),
    // array index: TR    | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,false,RowMajor,ColMajor,1>::run),
    // array index: ADJ   | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,Conj, RowMajor,ColMajor,1>::run),
    0,
    // array index: NOTR  | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,false,ColMajor,ColMajor,1>::run),
    // array index: TR    | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,false,RowMajor,ColMajor,1>::run),
    // array index: ADJ   | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,Conj, RowMajor,ColMajor,1>::run),
    0,
    // array index: NOTR  | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,false,ColMajor,ColMajor,1>::run),
    // array index: TR    | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,false,RowMajor,ColMajor,1>::run),
    // array index: ADJ   | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,Conj, RowMajor,ColMajor,1>::run),
    0
  };

  const Scalar* a = reinterpret_cast<const Scalar*>(pa);
  Scalar* b = reinterpret_cast<Scalar*>(pb);
  Scalar  alpha = *reinterpret_cast<const Scalar*>(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<std::max(1,(SIDE(*side)==LEFT)?*m:*n))                 info = 9;
  else if(*ldb<std::max(1,*m))                                        info = 11;
  if(info)
    return xerbla_(SCALAR_SUFFIX_UP"TRSM ",&info,6);

  if(*m==0 || *n==0)
    return 0;

  int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4);

  if(SIDE(*side)==LEFT)
  {
    internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*m,1,false);
    func[code](*m, *n, a, *lda, b, 1, *ldb, blocking);
  }
  else
  {
    internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*n,1,false);
    func[code](*n, *m, a, *lda, b, 1, *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)(const char *side, const char *uplo, const char *opa, const char *diag, const int *m, const int *n,
                          const RealScalar *palpha, const RealScalar *pa, const int *lda, RealScalar *pb, const 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, DenseIndex, const Scalar&, internal::level3_blocking<Scalar,Scalar>&);
  static const functype func[32] = {
    // array index: NOTR  | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          true, ColMajor,false,ColMajor,false,ColMajor,1>::run),
    // array index: TR    | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          true, RowMajor,false,ColMajor,false,ColMajor,1>::run),
    // array index: ADJ   | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          true, RowMajor,Conj, ColMajor,false,ColMajor,1>::run),
    0,
    // array index: NOTR  | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          false,ColMajor,false,ColMajor,false,ColMajor,1>::run),
    // array index: TR    | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          false,ColMajor,false,RowMajor,false,ColMajor,1>::run),
    // array index: ADJ   | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          false,ColMajor,false,RowMajor,Conj, ColMajor,1>::run),
    0,
    // array index: NOTR  | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          true, ColMajor,false,ColMajor,false,ColMajor,1>::run),
    // array index: TR    | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          true, RowMajor,false,ColMajor,false,ColMajor,1>::run),
    // array index: ADJ   | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          true, RowMajor,Conj, ColMajor,false,ColMajor,1>::run),
    0,
    // array index: NOTR  | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          false,ColMajor,false,ColMajor,false,ColMajor,1>::run),
    // array index: TR    | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          false,ColMajor,false,RowMajor,false,ColMajor,1>::run),
    // array index: ADJ   | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          false,ColMajor,false,RowMajor,Conj, ColMajor,1>::run),
    0,
    // array index: NOTR  | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor,1>::run),
    // array index: TR    | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor,1>::run),
    // array index: ADJ   | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor,1>::run),
    0,
    // array index: NOTR  | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor,1>::run),
    // array index: TR    | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor,1>::run),
    // array index: ADJ   | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor,1>::run),
    0,
    // array index: NOTR  | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor,1>::run),
    // array index: TR    | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor,1>::run),
    // array index: ADJ   | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor,1>::run),
    0,
    // array index: NOTR  | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor,1>::run),
    // array index: TR    | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor,1>::run),
    // array index: ADJ   | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor,1>::run),
    0
  };

  const Scalar* a = reinterpret_cast<const Scalar*>(pa);
  Scalar* b = reinterpret_cast<Scalar*>(pb);
  Scalar  alpha = *reinterpret_cast<const Scalar*>(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<std::max(1,(SIDE(*side)==LEFT)?*m:*n))                 info = 9;
  else if(*ldb<std::max(1,*m))                                        info = 11;
  if(info)
    return xerbla_(SCALAR_SUFFIX_UP"TRMM ",&info,6);

  int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4);

  if(*m==0 || *n==0)
    return 1;

  // FIXME find a way to avoid this copy
  Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp = matrix(b,*m,*n,*ldb);
  matrix(b,*m,*n,*ldb).setZero();

  if(SIDE(*side)==LEFT)
  {
    internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*m,1,false);
    func[code](*m, *n, *m, a, *lda, tmp.data(), tmp.outerStride(), b, 1, *ldb, alpha, blocking);
  }
  else
  {
    internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*n,1,false);
    func[code](*m, *n, *n, tmp.data(), tmp.outerStride(), a, *lda, b, 1, *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)(const char *side, const char *uplo, const int *m, const int *n, const RealScalar *palpha,
                          const RealScalar *pa, const int *lda, const RealScalar *pb, const int *ldb, const RealScalar *pbeta, RealScalar *pc, const int *ldc)
{
//   std::cerr << "in symm " << *side << " " << *uplo << " " << *m << "x" << *n << " lda:" << *lda << " ldb:" << *ldb << " ldc:" << *ldc << " alpha:" << *palpha << " beta:" << *pbeta << "\n";
  const Scalar* a = reinterpret_cast<const Scalar*>(pa);
  const Scalar* b = reinterpret_cast<const Scalar*>(pb);
  Scalar* c = reinterpret_cast<Scalar*>(pc);
  Scalar alpha = *reinterpret_cast<const Scalar*>(palpha);
  Scalar beta  = *reinterpret_cast<const Scalar*>(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<std::max(1,(SIDE(*side)==LEFT)?*m:*n))                 info = 7;
  else if(*ldb<std::max(1,*m))                                        info = 9;
  else if(*ldc<std::max(1,*m))                                        info = 12;
  if(info)
    return xerbla_(SCALAR_SUFFIX_UP"SYMM ",&info,6);

  if(beta!=Scalar(1))
  {
    if(beta==Scalar(0)) matrix(c, *m, *n, *ldc).setZero();
    else                matrix(c, *m, *n, *ldc) *= beta;
  }

  if(*m==0 || *n==0)
  {
    return 1;
  }

  int size = (SIDE(*side)==LEFT) ? (*m) : (*n);
  #if ISCOMPLEX
  // FIXME add support for symmetric complex matrix
  Matrix<Scalar,Dynamic,Dynamic,ColMajor> matA(size,size);
  if(UPLO(*uplo)==UP)
  {
    matA.triangularView<Upper>() = matrix(a,size,size,*lda);
    matA.triangularView<Lower>() = matrix(a,size,size,*lda).transpose();
  }
  else if(UPLO(*uplo)==LO)
  {
    matA.triangularView<Lower>() = matrix(a,size,size,*lda);
    matA.triangularView<Upper>() = 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
  internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic> blocking(*m,*n,size,1,false);

  if(SIDE(*side)==LEFT)
    if(UPLO(*uplo)==UP)       internal::product_selfadjoint_matrix<Scalar, DenseIndex, RowMajor,true,false, ColMajor,false,false, ColMajor,1>::run(*m, *n, a, *lda, b, *ldb, c, 1, *ldc, alpha, blocking);
    else if(UPLO(*uplo)==LO)  internal::product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor,true,false, ColMajor,false,false, ColMajor,1>::run(*m, *n, a, *lda, b, *ldb, c, 1, *ldc, alpha, blocking);
    else                      return 0;
  else if(SIDE(*side)==RIGHT)
    if(UPLO(*uplo)==UP)       internal::product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor,false,false, RowMajor,true,false, ColMajor,1>::run(*m, *n, b, *ldb, a, *lda, c, 1, *ldc, alpha, blocking);
    else if(UPLO(*uplo)==LO)  internal::product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor,false,false, ColMajor,true,false, ColMajor,1>::run(*m, *n, b, *ldb, a, *lda, c, 1, *ldc, alpha, blocking);
    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)(const char *uplo, const char *op, const int *n, const int *k,
                          const RealScalar *palpha, const RealScalar *pa, const int *lda, const RealScalar *pbeta, RealScalar *pc, const 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, DenseIndex, const Scalar&, internal::level3_blocking<Scalar,Scalar>&);
  static const functype func[8] = {
    // array index: NOTR  | (UP << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,ColMajor,Conj, 1, Upper>::run),
    // array index: TR    | (UP << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,ColMajor,Conj, 1, Upper>::run),
    // array index: ADJ   | (UP << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,ColMajor,false,1, Upper>::run),
    0,
    // array index: NOTR  | (LO << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,ColMajor,Conj, 1, Lower>::run),
    // array index: TR    | (LO << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,ColMajor,Conj, 1, Lower>::run),
    // array index: ADJ   | (LO << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,ColMajor,false,1, Lower>::run),
    0
  };
  #endif

  const Scalar* a = reinterpret_cast<const Scalar*>(pa);
  Scalar* c = reinterpret_cast<Scalar*>(pc);
  Scalar alpha = *reinterpret_cast<const Scalar*>(palpha);
  Scalar beta  = *reinterpret_cast<const Scalar*>(pbeta);

  int info = 0;
  if(UPLO(*uplo)==INVALID)                                            info = 1;
  else if(OP(*op)==INVALID || (ISCOMPLEX && OP(*op)==ADJ) )           info = 2;
  else if(*n<0)                                                       info = 3;
  else if(*k<0)                                                       info = 4;
  else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 7;
  else if(*ldc<std::max(1,*n))                                        info = 10;
  if(info)
    return xerbla_(SCALAR_SUFFIX_UP"SYRK ",&info,6);

  if(beta!=Scalar(1))
  {
    if(UPLO(*uplo)==UP)
      if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
      else                matrix(c, *n, *n, *ldc).triangularView<Upper>() *= beta;
    else
      if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
      else                matrix(c, *n, *n, *ldc).triangularView<Lower>() *= beta;
  }

  if(*n==0 || *k==0)
    return 0;

  #if ISCOMPLEX
  // FIXME add support for symmetric complex matrix
  if(UPLO(*uplo)==UP)
  {
    if(OP(*op)==NOTR)
      matrix(c, *n, *n, *ldc).triangularView<Upper>() += alpha * matrix(a,*n,*k,*lda) * matrix(a,*n,*k,*lda).transpose();
    else
      matrix(c, *n, *n, *ldc).triangularView<Upper>() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda);
  }
  else
  {
    if(OP(*op)==NOTR)
      matrix(c, *n, *n, *ldc).triangularView<Lower>() += alpha * matrix(a,*n,*k,*lda) * matrix(a,*n,*k,*lda).transpose();
    else
      matrix(c, *n, *n, *ldc).triangularView<Lower>() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda);
  }
  #else
  internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic> blocking(*n,*n,*k,1,false);

  int code = OP(*op) | (UPLO(*uplo) << 2);
  func[code](*n, *k, a, *lda, a, *lda, c, 1, *ldc, alpha, blocking);
  #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)(const char *uplo, const char *op, const int *n, const int *k, const RealScalar *palpha,
                           const RealScalar *pa, const int *lda, const RealScalar *pb, const int *ldb, const RealScalar *pbeta, RealScalar *pc, const int *ldc)
{
  const Scalar* a = reinterpret_cast<const Scalar*>(pa);
  const Scalar* b = reinterpret_cast<const Scalar*>(pb);
  Scalar* c = reinterpret_cast<Scalar*>(pc);
  Scalar alpha = *reinterpret_cast<const Scalar*>(palpha);
  Scalar beta  = *reinterpret_cast<const Scalar*>(pbeta);

//   std::cerr << "in syr2k " << *uplo << " " << *op << " " << *n << " " << *k << " " << alpha << " " << *lda << " " << *ldb << " " << beta << " " << *ldc << "\n";

  int info = 0;
  if(UPLO(*uplo)==INVALID)                                            info = 1;
  else if(OP(*op)==INVALID || (ISCOMPLEX && OP(*op)==ADJ) )           info = 2;
  else if(*n<0)                                                       info = 3;
  else if(*k<0)                                                       info = 4;
  else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 7;
  else if(*ldb<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 9;
  else if(*ldc<std::max(1,*n))                                        info = 12;
  if(info)
    return xerbla_(SCALAR_SUFFIX_UP"SYR2K",&info,6);

  if(beta!=Scalar(1))
  {
    if(UPLO(*uplo)==UP)
      if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
      else                matrix(c, *n, *n, *ldc).triangularView<Upper>() *= beta;
    else
      if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
      else                matrix(c, *n, *n, *ldc).triangularView<Lower>() *= beta;
  }

  if(*k==0)
    return 1;

  if(OP(*op)==NOTR)
  {
    if(UPLO(*uplo)==UP)
    {
      matrix(c, *n, *n, *ldc).triangularView<Upper>()
        += 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<Lower>()
        += 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<Upper>()
        += 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<Lower>()
        += 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)(const char *side, const char *uplo, const int *m, const int *n, const RealScalar *palpha,
                          const RealScalar *pa, const int *lda, const RealScalar *pb, const int *ldb, const RealScalar *pbeta, RealScalar *pc, const int *ldc)
{
  const Scalar* a = reinterpret_cast<const Scalar*>(pa);
  const Scalar* b = reinterpret_cast<const Scalar*>(pb);
  Scalar* c = reinterpret_cast<Scalar*>(pc);
  Scalar alpha = *reinterpret_cast<const Scalar*>(palpha);
  Scalar beta  = *reinterpret_cast<const Scalar*>(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<std::max(1,(SIDE(*side)==LEFT)?*m:*n))                 info = 7;
  else if(*ldb<std::max(1,*m))                                        info = 9;
  else if(*ldc<std::max(1,*m))                                        info = 12;
  if(info)
    return xerbla_(SCALAR_SUFFIX_UP"HEMM ",&info,6);

  if(beta==Scalar(0))       matrix(c, *m, *n, *ldc).setZero();
  else if(beta!=Scalar(1))  matrix(c, *m, *n, *ldc) *= beta;

  if(*m==0 || *n==0)
  {
    return 1;
  }

  int size = (SIDE(*side)==LEFT) ? (*m) : (*n);
  internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic> blocking(*m,*n,size,1,false);

  if(SIDE(*side)==LEFT)
  {
    if(UPLO(*uplo)==UP)       internal::product_selfadjoint_matrix<Scalar,DenseIndex,RowMajor,true,Conj,  ColMajor,false,false, ColMajor, 1>
                                ::run(*m, *n, a, *lda, b, *ldb, c, 1, *ldc, alpha, blocking);
    else if(UPLO(*uplo)==LO)  internal::product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,true,false, ColMajor,false,false, ColMajor,1>
                                ::run(*m, *n, a, *lda, b, *ldb, c, 1, *ldc, alpha, blocking);
    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<Upper>();/*internal::product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,false,false, RowMajor,true,Conj,  ColMajor, 1>
                                ::run(*m, *n, b, *ldb, a, *lda, c, 1, *ldc, alpha, blocking);*/
    else if(UPLO(*uplo)==LO)  internal::product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,false,false, ColMajor,true,false, ColMajor,1>
                                ::run(*m, *n, b, *ldb, a, *lda, c, 1, *ldc, alpha, blocking);
    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)(const char *uplo, const char *op, const int *n, const int *k,
                          const RealScalar *palpha, const RealScalar *pa, const int *lda, const RealScalar *pbeta, RealScalar *pc, const int *ldc)
{
//   std::cerr << "in herk " << *uplo << " " << *op << " " << *n << " " << *k << " " << *palpha << " " << *lda << " " << *pbeta << " " << *ldc << "\n";

  typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, DenseIndex, const Scalar&, internal::level3_blocking<Scalar,Scalar>&);
  static const functype func[8] = {
    // array index: NOTR  | (UP << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,1,Upper>::run),
    0,
    // array index: ADJ   | (UP << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,1,Upper>::run),
    0,
    // array index: NOTR  | (LO << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,1,Lower>::run),
    0,
    // array index: ADJ   | (LO << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,1,Lower>::run),
    0
  };

  const Scalar* a = reinterpret_cast<const Scalar*>(pa);
  Scalar* c = reinterpret_cast<Scalar*>(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<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 7;
  else if(*ldc<std::max(1,*n))                                        info = 10;
  if(info)
    return xerbla_(SCALAR_SUFFIX_UP"HERK ",&info,6);

  int code = OP(*op) | (UPLO(*uplo) << 2);

  if(beta!=RealScalar(1))
  {
    if(UPLO(*uplo)==UP)
      if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
      else                matrix(c, *n, *n, *ldc).triangularView<StrictlyUpper>() *= beta;
    else
      if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
      else                matrix(c, *n, *n, *ldc).triangularView<StrictlyLower>() *= 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))
  {
    internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic> blocking(*n,*n,*k,1,false);
    func[code](*n, *k, a, *lda, a, *lda, c, 1, *ldc, alpha, blocking);
    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)(const char *uplo, const char *op, const int *n, const int *k,
                           const RealScalar *palpha, const RealScalar *pa, const int *lda, const RealScalar *pb, const int *ldb, const RealScalar *pbeta, RealScalar *pc, const int *ldc)
{
  const Scalar* a = reinterpret_cast<const Scalar*>(pa);
  const Scalar* b = reinterpret_cast<const Scalar*>(pb);
  Scalar* c = reinterpret_cast<Scalar*>(pc);
  Scalar alpha = *reinterpret_cast<const Scalar*>(palpha);
  RealScalar beta  = *pbeta;

//   std::cerr << "in her2k " << *uplo << " " << *op << " " << *n << " " << *k << " " << alpha << " " << *lda << " " << *ldb << " " << 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<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 7;
  else if(*ldb<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 9;
  else if(*ldc<std::max(1,*n))                                        info = 12;
  if(info)
    return xerbla_(SCALAR_SUFFIX_UP"HER2K",&info,6);

  if(beta!=RealScalar(1))
  {
    if(UPLO(*uplo)==UP)
      if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
      else                matrix(c, *n, *n, *ldc).triangularView<StrictlyUpper>() *= beta;
    else
      if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
      else                matrix(c, *n, *n, *ldc).triangularView<StrictlyLower>() *= 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<Upper>()
        +=            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<Lower>()
        += 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<Upper>()
        +=             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<Lower>()
        +=             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

quality100%

¤ Dauer der Verarbeitung: 0.22 Sekunden (vorverarbeitet) ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.